MODULE exchtestmod
  USE par_kind,     ONLY: wp
  USE mapcomm_mod,  ONLY: cyclic_bc, trimmed, nidx, eidx, sidx, widx
  PRIVATE

  LOGICAL, PARAMETER :: stop_after_testing  = .FALSE., & ! STOP the nemo run 
                                                         ! once these tests are complete
#if defined key_mpp_rkpart
                        compare_exch_methods= .FALSE., & ! Do both packed and 
                        do_integer_tests    = .FALSE., & ! indiv. exchanges and compare output
                        do_real_tests       = .TRUE. , &
                        do_integer_timings  = .FALSE., &
                        do_real_timings     = .TRUE.,  &
                        use_exch_list       = .FALSE.   ! Whether to use the halo
                                                        ! packing API for the tests
                                                        ! - NOT currently working!
#else
                        compare_exch_methods= .FALSE., & 
                      ! WARNING: test code not supported for non rkpart build!
                        do_integer_tests    = .FALSE., & 
                        do_real_tests       = .FALSE., &
                        do_integer_timings  = .FALSE., &
                        do_real_timings     = .FALSE., &
                        use_exch_list       = .FALSE.   ! Whether to use the halo
                                                        ! packing API for the tests
                                                        ! - NOT currently working!
#endif

  LOGICAL :: test_failed  ! Updated if any test fails so we can
                          ! stop the job once testing complete

  ! We need to test exchanges of {REAL, INTEGER} {2D, 3D} arrays on
  ! {U, V, W, Z, T} grids
!FTRANS r3d   :I :I :z
!FTRANS r3d_2 :I :I :z
!FTRANS r3d_3 :I :I :z
  REAL (kind=wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: r3d, r3d_2, r3d_3
  REAL (kind=wp), DIMENSION(:,:)  , ALLOCATABLE, TARGET :: r2d, r2d_2
  REAL (kind=wp), DIMENSION(:,:)  , ALLOCATABLE         :: r2dxz
!FTRANS i3d   :I :I :z
!FTRANS i3d_2 :I :I :z
  INTEGER, DIMENSION(:,:,:), ALLOCATABLE, TARGET :: i3d, i3d_2
  INTEGER, DIMENSION(:,:),   ALLOCATABLE, TARGET :: i2d, i2d_2

  ! Last ocean level above ocean floor
  INTEGER, DIMENSION(:,:),   POINTER             :: pmaxdepth

  ! Unit to use for outputting log of results
  INTEGER, PARAMETER :: LOG_UNIT = 1002

  ! Identifiers for types of test
  INTEGER, PARAMETER :: PE_ID_TEST       = 0, &
                        GLOBAL_LOCN_TEST = 1

  ! No. of times to repeat a message exchange while doing timing
  INTEGER, PARAMETER :: loops_to_time = 5000

  INTERFACE compare_arrays
    MODULE PROCEDURE compare_arrays2d, compare_arrays3d
  END INTERFACE compare_arrays

  PUBLIC mpp_test_comms ! Only the top-level driver routine is public

CONTAINS

  !====================================================================

  SUBROUTINE mpp_test_comms(depth, lmaxdepth)
    USE par_oce,      ONLY: jpi, jpj, jpk, jpkorig, jpreci
    USE par_kind,     ONLY: wp
    USE mapcomm_mod,  ONLY: Iminus, Iplus, Jminus, Jplus, NONE
    USE lib_mpp,      ONLY: ctl_stop, ctl_warn
    USE dom_oce,      ONLY: narea
    USE exchmod,      ONLY: add_exch, bound_exch_list, bound_exch
    USE timing,       ONLY: timing_finalize
    USE profile
#if defined key_mpp_mpi
    USE lib_mpp,      ONLY: mpi_comm_opa
    USE mpi
#endif
    IMPLICIT none
    ! Routine arguments
    INTEGER, DIMENSION(:,:), INTENT(in) :: depth    ! Mask (1 for ocean, 0 for land)
    INTEGER, DIMENSION(:,:), TARGET     :: lmaxdepth ! Last level above ocean floor
    ! Local vars
    INTEGER            :: ierr
    CHARACTER(len=256) :: name

#if ! defined key_mpp_mpi
    CALL ctl_warn('mpp_test_comms: not built with MPI so nothing to do!')
    RETURN
#endif

    CALL prof_tracing_on()

    ! Initialise the flag that will store whether we have any fails
    test_failed = .FALSE.

    ALLOCATE(r3d(jpi,jpj,jpkorig), r3d_2(jpi,jpj,jpkorig), r3d_3(jpi,jpj,jpkorig), &
             r2d(jpi,jpj), r2d_2(jpi,jpj), Stat=ierr)

    ! Set module member variable to point to max-depth data so we can access it
    ! when checking results of halo swaps.
    pmaxdepth => lmaxdepth

    IF(ierr .ne. 0)THEN
       WRITE (*,*) "Failed to allocate memory in mpp_test_comms - no tests will be performed!"
       RETURN
    END IF

    WRITE(name, FMT="('testing.',I4.4,'.log')") narea-1
    OPEN(UNIT=LOG_UNIT, FILE=TRIM(name), &
         STATUS='REPLACE', ACTION='WRITE', IOSTAT=ierr)

    IF(ierr .ne. 0)THEN
       WRITE (*,*) "Failed to open testing.log for logging - no tests will be performed!"
       RETURN
    END IF

    ! Test exchanges of real arrays ------------------------------------------

    IF(do_real_tests)THEN

       IF(narea == 1) WRITE (*,*) '1. Testing exchanges of real arrays...'

       ! 1. Test halo exchanges for a 3D REAL array at 'W' point...

       name = '3D REAL array at W point'
       CALL exch_test(name, 'W', GLOBAL_LOCN_TEST, depth, r3d1=r3d, isgn=-1)

       name = '3D REAL array at W point, lfill reversed'
       CALL exch_test(name, 'W', GLOBAL_LOCN_TEST, depth, r3d1=r3d, &
                      lfill=.TRUE.)

       IF(narea == 1) WRITE (*,*) 'Test 1 done.'

       ! 2. Test halo exchanges for a 3D REAL array at 'Z' point...

       name = '3D REAL array at Z point'
       CALL exch_test(name, 'Z', GLOBAL_LOCN_TEST, depth, &
                       r3d1=r3d, isgn=-999, lfill=.TRUE.)

       IF(narea == 1) WRITE (*,*) 'Test 2 done.'

       ! 3. Test halo exchanges for a 2D REAL array at 'Z' point...

       name = '2D REAL array at Z point'
       CALL exch_test(name, 'Z', GLOBAL_LOCN_TEST, depth, r2d1=r2d, &
                      isgn=-999, lfill=.TRUE.)

       IF(narea == 1) WRITE (*,*) 'Test 3 done.'

       ! 4. Test halo exchanges for a 2D REAL array at 'T' point...

       name = '2D REAL array at T point'
       CALL exch_test(name, 'T', GLOBAL_LOCN_TEST, depth, r2d1=r2d, &
                      isgn=1, lfill=.FALSE.)

       IF(narea == 1) WRITE (*,*) 'Test 4 done.'


       ! 5. Test halo exchanges for a 3D REAL array at 'T' point...

       name = '3D REAL array at T point'
       CALL exch_test(name, 'T', GLOBAL_LOCN_TEST, depth, r3d1=r3d, &
                      isgn=-999, lfill=.TRUE.)

       name = '3D REAL array at T point, lfill now .FALSE.'
       CALL exch_test(name, 'T', GLOBAL_LOCN_TEST, depth, r3d1=r3d)

       IF(narea == 1) WRITE (*,*) 'Test 5 done.'

       ! 6. Test halo exchanges for a 3D REAL array at 'V' point...

       name = '3D REAL array at V point'
       CALL exch_test(name, 'V', PE_ID_TEST, depth, r3d1=r3d, &
            isgn=-999, lfill=.TRUE.)

       IF(narea == 1) WRITE (*,*) 'Test 6 done.'

       ! 7. Test halo exchanges for _two_ 2D REAL arrays at 'V' point...

       name = 'Two 2D REAL arrays at V point'
       CALL exch_test(name,'V', PE_ID_TEST, depth, r2d1=r2d, r2d2=r2d_2, &
            isgn=-999, lfill=.TRUE.)

       IF(narea == 1) WRITE (*,*) 'Test 7 done.'

       ! 8. Test halo exchanges for _two_ 3D REAL arrays at 'V' point...

       name = 'Two 3D REAL arrays at V point'
       CALL exch_test(name,'V', PE_ID_TEST, depth, r3d1=r3d, r3d2=r3d_2, &
            isgn=-999, lfill=.TRUE.)

       IF(narea == 1) WRITE (*,*) 'Test 8 done.'

       ! 9. Test halo exchanges for _three_ 3D REAL arrays at 'T' point...

       name = 'Three 3D REAL arrays at T point'
       CALL exch_test(name,'T', PE_ID_TEST, depth,      &
                      r3d1=r3d, r3d2=r3d_2, r3d3=r3d_3, &
                      isgn=-999, lfill=.TRUE.)

       IF(narea == 1) WRITE (*,*) 'Test 9 done.'

    END IF

    ! End of real-array tests, now for timing ------------------------------

    IF(do_real_timings)THEN

       IF(narea == 1)THEN
          WRITE (*,*) '2. Measuring performance for real arrays...'
       END IF

       CALL exch_timing('Z', r2d1=r2d,             &
                        isgn=-999, lfill=.TRUE.)
       CALL exch_timing('T', r2d1=r2d, r2d2=r2d_2, &
                        isgn=-999, lfill=.TRUE.)
       CALL exch_timing('T', r2d1=r2d)
       CALL exch_timing('U', r2d1=r2d, isgn=-1)
       CALL exch_timing('T', r2d1=r2d, r2d2=r2d_2)

       CALL exch_timing('T', r3d1=r3d)
       CALL exch_timing('Z', r3d1=r3d, r2d1=r2d,   &
                        isgn=-999, lfill=.TRUE.)
       CALL exch_timing('U', r3d1=r3d, r3d2=r3d_2, &
                        isgn=-1, lfill=.TRUE.)
       CALL exch_timing('U', r3d1=r3d, r3d2=r3d_2, isgn=-1)
       CALL exch_timing('T', r3d1=r3d, r3d2=r3d_2, &
                        isgn=-1, lfill=.TRUE.)
       CALL exch_timing('T', r3d1=r3d, r3d2=r3d_2)

       IF(narea == 1)THEN
          WRITE (*,*) '...performance tests for real array exchanges complete.'
       END IF

    END IF

    ! Free-up memory used for real-array tests

    DEALLOCATE(r3d, r3d_2, r3d_3, r2d)
    ! and allocate that for integer-array tests
    ALLOCATE(i3d(jpi,jpj,jpkorig), i3d_2(jpi,jpj,jpkorig), &
             i2d(jpi,jpj), i2d_2(jpi,jpj), Stat=ierr)

    IF(ierr .ne. 0)THEN
       WRITE (*,*) "Failed to allocate memory in mpp_test_comms - no INTEGER tests will be performed!"
       RETURN
    END IF

    ! Test exchanges of integer arrays ------------------------------------------

    IF(do_integer_tests)THEN

       IF(narea == 1)THEN
          WRITE (*,*) '3. Testing exchanges of integer arrays...'
       END IF

       ! 8. Test halo exchanges for a 3D INTEGER array at 'W' point...

       name = 'One 3D integer array at W point'
       CALL exch_test(name,'W', PE_ID_TEST, depth, i3d1=i3d, &
            isgn=-999, lfill=.TRUE.)

       ! 9. Test halo exchanges for a 3D INTEGER array at 'Z' point...

       name = 'One 3D integer array at Z point'
       CALL exch_test(name,'Z', PE_ID_TEST, depth, i3d1=i3d, &
            isgn=-999, lfill=.TRUE.)

       ! 10. Test halo exchanges for a 2D INTEGER array at 'Z' point...

       name = 'One 2D integer array at Z point'
       CALL exch_test(name,'Z', PE_ID_TEST, depth, i2d1=i2d, &
            isgn=-999, lfill=.TRUE.)

       ! 11. Test halo exchanges for _two_ 3D INTEGER arrays at 'Z' point...

       name = 'Two 3D integer arrays at Z point'
       CALL exch_test(name,'Z', PE_ID_TEST, depth, i3d1=i3d, i3d2=i3d_2, &
            isgn=-999, lfill=.TRUE.)


       IF(narea == 1)THEN
          WRITE (*,*) '...integer-array tests complete'
       END IF

    END IF

    ! End of tests, now for timing ------------------------------------

    IF(do_integer_timings)THEN

       IF(narea == 1)THEN
          WRITE (*,*) '4. Measuring performance for integer arrays...'
       END IF

       CALL exch_timing('Z', i2d1=i2d, &
            isgn=-999, lfill=.TRUE.)

       CALL exch_timing('Z', i2d1=i2d, i2d2=i2d_2, &
            isgn=-999, lfill=.TRUE.)

       CALL exch_timing('Z', i3d1=i3d, &
            isgn=-999, lfill=.TRUE.)

       CALL exch_timing('Z', i2d1 = i2d, i3d1=i3d, &
            isgn=-999, lfill=.TRUE.)

       CALL exch_timing('Z', i3d1=i3d, i3d2=i3d_2, &
            isgn=-999, lfill=.TRUE.)

       CALL exch_timing('Z', i2d1 = i2d, i3d1=i3d, i3d2=i3d_2, &
            isgn=-999, lfill=.TRUE.)

       IF(narea == 1)THEN
          WRITE (*,*) '...integer-array performance tests complete.'
       END IF

    END IF

    ! End of timing of integer exchanges -------------------------------

    CLOSE(UNIT=LOG_UNIT)

    DEALLOCATE(i3d, i3d_2, i2d, i2d_2)

#if defined key_mpp_mpi
    ! Check for success or otherwise of tests on all PEs
    CALL mpi_allreduce(MPI_IN_PLACE, test_failed, 1, MPI_LOGICAL, MPI_LOR, &
                       mpi_comm_opa, ierr )
#endif

    IF(stop_after_testing .OR. test_failed )THEN
       IF(narea == 1)THEN
          IF(test_failed)THEN
             WRITE (*,FMT="('Stopping due to error in msg. exchange tests!')")
          ELSE
             WRITE (*,FMT="('Stopping now that comms tests are complete!')")
          END IF
       END IF
       ! Generate a timing report
       CALL timing_finalize()
       ! Dirty way of causing NEMO to stop immediately
       CALL ctl_stop('STOP', 'Stopping now that comms tests are complete')
    END IF

  END SUBROUTINE mpp_test_comms

  !======================================================================

  SUBROUTINE exch_test(descr, gridType, testType, depth, r2d1, r2d2, &
                       r3d1, r3d2, r3d3, i2d1, i3d1, i3d2, isgn, lfill)
    USE par_kind,     ONLY: wp
    USE par_oce,      ONLY: jpreci
    USE mapcomm_mod,  ONLY: Iminus, Iplus, Jminus, Jplus, NONE
    USE exchmod,      ONLY: add_exch, bound_exch_list, bound_exch
    USE lbclnk,       ONLY: lbc_lnk
    USE lib_mpp,      ONLY: ctl_warn
    USE dom_oce,      ONLY: narea
    IMPLICIT none
!FTRANS r3d1 :I :I :z
!FTRANS r3d2 :I :I :z
!FTRANS r3d3 :I :I :z
!FTRANS i3d1 :I :I :z
!FTRANS i3d2 :I :I :z
!FTRANS r3dcopy :I :I :z
    ! Arguments
    CHARACTER(LEN=256), INTENT(in) :: descr    ! Description of the test being done
    CHARACTER(LEN=1),   INTENT(in) :: gridType ! Grid on which test being done
    INTEGER,            INTENT(in) :: testType ! Type of test to do
    INTEGER, DIMENSION(:,:),          INTENT(in)              :: depth ! Global land mask
    LOGICAL, INTENT(in),                             OPTIONAL :: lfill
    REAL (kind=wp), DIMENSION(:,:),   INTENT(inout), OPTIONAL :: r2d1,r2d2
    REAL (kind=wp), DIMENSION(:,:,:), INTENT(inout), OPTIONAL :: r3d1,r3d2,r3d3
    INTEGER, DIMENSION(:,:,:),        INTENT(inout), OPTIONAL :: i3d1,i3d2
    INTEGER, DIMENSION(:,:),          INTENT(inout), OPTIONAL :: i2d1
    INTEGER,                          INTENT(in),    OPTIONAL :: isgn
    ! Local vars
    LOGICAL       :: hit_error, test_pass
    REAL (kind=wp), DIMENSION(:,:),   ALLOCATABLE :: r2dcopy
    REAL (kind=wp), DIMENSION(:,:,:), ALLOCATABLE :: r3dcopy
    INTEGER       :: ierr
    INTEGER       :: isize_x, isize_y, isize_z
    REAL (kind=wp):: psgn
    !!-----------------------------------------------------------------------

#if ! defined key_mpp_rkpart
    CALL ctl_warn('exch_test: halo exchange testing not supported for build without key_mpp_rkpart defined')
    RETURN
#endif

    ! Initialise arrays being exchanged
    ! A correct exchange process (but without north-fold) won't change
    ! these values. If lfill is set to .TRUE. then no north-fold exchange
    ! is performed.
    CALL init_test_arrays(testType, &
                          r2d1, r2d2, r3d1, r3d2, r3d3, i2d1, i3d1, i3d2)
       
    IF(PRESENT(r3d1))THEN

#if defined key_mpp_rkpart
       IF(use_exch_list)THEN
          CALL add_exch(jpreci, gridType, Iminus,Jminus,Iplus,Jplus, r3d=r3d1, &
                        isgn=isgn, lfill=lfill)
       ELSE
          CALL bound_exch(r3d1, jpreci, jpreci,      &
                          Iplus, Iminus, Jplus, Jminus, &
                          gridType, isgn=isgn, lfill=lfill)
       END IF
#else
       ! psgn is NOT an optional argument to lbc_lnk
       psgn = 1.0_wp
       IF(PRESENT(isgn))psgn = REAL(isgn)
       CALL lbc_lnk(r3d1, gridType, psgn)
#endif
    END IF

#if defined key_mpp_rkpart
    IF(PRESENT(r3d2))THEN
       IF(use_exch_list)THEN
          CALL add_exch(jpreci, gridType, Iminus,Jminus,Iplus,Jplus, r3d=r3d2, &
                        isgn=isgn, lfill=lfill)
       ELSE
          CALL bound_exch(r3d2, jpreci, jpreci,      &
                          Iplus, Iminus, Jplus, Jminus, &
                          gridType, isgn=isgn, lfill=lfill)
       END IF
    END IF
#endif

    IF(PRESENT(r3d3))THEN

#if defined key_mpp_rkpart
       IF(use_exch_list)THEN
          CALL add_exch(jpreci, gridType, Iminus,Jminus,Iplus,Jplus, r3d=r3d3, &
                        isgn=isgn, lfill=lfill)
       ELSE
          CALL bound_exch(r3d3, jpreci, jpreci,      &
                          Iplus, Iminus, Jplus, Jminus, &
                          gridType, isgn=isgn, lfill=lfill)
       END IF
#endif

    END IF

    IF(PRESENT(r2d1))THEN

#if defined key_mpp_rkpart
       IF(use_exch_list)THEN
          CALL add_exch(jpreci, gridType, Iminus,Jminus,Iplus,Jplus, r2d=r2d1, &
                        isgn=isgn, lfill=lfill)
       ELSE
          CALL bound_exch(r2d1, jpreci, jpreci,      &
                          Iplus, Iminus, Jplus, Jminus, &
                          gridType, isgn=isgn, lfill=lfill)
       END IF
#endif
    END IF

    IF(PRESENT(r2d2))THEN

#if defined key_mpp_rkpart
       IF(use_exch_list)THEN
          CALL add_exch(jpreci, gridType, Iminus,Jminus,Iplus,Jplus, r2d=r2d2, &
                        isgn=isgn, lfill=lfill)
       ELSE
          CALL bound_exch(r2d2, jpreci, jpreci,      &
                          Iplus, Iminus, Jplus, Jminus, &
                          gridType, isgn=isgn, lfill=lfill)
       END IF
#endif
    END IF

    IF( PRESENT(i2d1) )THEN

#if defined key_mpp_rkpart
       IF(use_exch_list)THEN
          CALL add_exch(jpreci, gridType, Iminus,Jminus,Iplus,Jplus, i2d=i2d1, &
                        isgn=isgn, lfill=lfill)
        ELSE
          CALL bound_exch(i2d1, jpreci, jpreci,      &
                          Iplus, Iminus, Jplus, Jminus, &
                          gridType, isgn=isgn, lfill=lfill)
       END IF
#endif
   END IF

    IF( PRESENT(i3d1) )THEN

#if defined key_mpp_rkpart
       IF(use_exch_list)THEN
       CALL add_exch(jpreci, gridType, Iminus,Jminus,Iplus,Jplus, i3d=i3d1, &
                     isgn=isgn, lfill=lfill)
       ELSE
          CALL bound_exch(i3d1, jpreci, jpreci,      &
                          Iplus, Iminus, Jplus, Jminus, &
                          gridType, isgn=isgn, lfill=lfill)
       END IF
#endif
    END IF

    IF( PRESENT(i3d2) )THEN

#if defined key_mpp_rkpart
       IF(use_exch_list)THEN
       CALL add_exch(jpreci, gridType, Iminus,Jminus,Iplus,Jplus, i3d=i3d2, &
                     isgn=isgn, lfill=lfill)
       ELSE
          CALL bound_exch(i3d2, jpreci, jpreci,      &
                          Iplus, Iminus, Jplus, Jminus, &
                          gridType, isgn=isgn, lfill=lfill)
       END IF
#endif
    END IF
   
#if defined key_mpp_rkpart
    IF(use_exch_list)THEN
       ! Finally, exchange halos for the arrays that we've added to the list
       CALL bound_exch_list()
    END IF
#endif

    ! In this section we test to see that the results of using the halo-packing
    ! API are the same as doing individual halo swaps.
#if defined key_mpp_rkpart
    IF(compare_exch_methods .AND. use_exch_list)THEN

       WRITE (LOG_UNIT, FMT="('Comparing array exch. methods...')")

       IF(PRESENT(r3d1))THEN

#if defined key_z_first
          ! Store the results of the SIZE calls in variables otherwise ftrans falls
          ! over when dealing with the ALLOCATE call.
          isize_z = SIZE(r3d1,1)
          isize_y = SIZE(r3d1,3)
          isize_x =  SIZE(r3d1,2)
          ALLOCATE(r3dcopy(isize_x, isize_y, isize_z), Stat=ierr)
#else
          ALLOCATE(r3dcopy(SIZE(r3d1,1),SIZE(r3d1,2),SIZE(r3d1,3)), Stat=ierr) 
#endif
          IF(ierr == 0)THEN
             CALL init_test_arrays(testType, r3d1=r3dcopy)

             CALL bound_exch(r3dcopy, jpreci, jpreci,      &
                             Iplus, Iminus, Jplus, Jminus, &
                             gridType, isgn=isgn, lfill=lfill)

             CALL compare_arrays(r3d1, r3dcopy)
          ELSE
             WRITE (LOG_UNIT, FMT="('ERROR: failed to allocate r3dcopy to compare array exch. methods.')")
          END IF
       END IF

       IF(PRESENT(r2d1))THEN
          ALLOCATE(r2dcopy(SIZE(r2d1,1),SIZE(r2d1,2)), Stat=ierr) 

          IF(ierr == 0)THEN

             CALL init_test_arrays(testType, r2d1=r2dcopy)

             CALL bound_exch(r2dcopy, jpreci, jpreci,   &
                             Iminus,Jminus,Iplus,Jplus, &
                             gridType, isgn=isgn, lfill=lfill)

             CALL compare_arrays(r2d1, r2dcopy)
          ELSE
             WRITE (LOG_UNIT, FMT="('ERROR: failed to allocate r2dcopy to compare array exch. methods.')")
          END IF
       END IF
    END IF
#endif

    ! Use test_pass as a temporary, logical variable here
    test_pass = .FALSE.
    IF(PRESENT(lfill))test_pass = lfill

    SELECT CASE(testType)

       CASE (PE_ID_TEST)
          WRITE(LOG_UNIT, FMT="('Testing results of exchanging PE IDs')")
       CASE(GLOBAL_LOCN_TEST)
          WRITE(LOG_UNIT, FMT="('Testing results of exchanging global coords')")
       CASE DEFAULT
          WRITE(LOG_UNIT, FMT="('ERROR: unrecognised test type in exch_test()')")
          test_pass = .FALSE.
          RETURN
    END SELECT
    WRITE(LOG_UNIT, FMT="((A),' with lfill=',L3,':')") TRIM(descr), test_pass

    test_pass = .TRUE.

    IF( PRESENT(r3d1) )THEN
       CALL array_check(gridType, testType, depth, r3d=r3d1, lfill=lfill, &
                        stat=hit_error)
       IF(test_pass)test_pass = .NOT. hit_error
    END IF
    IF( PRESENT(r3d2) )THEN
       CALL array_check(gridType, testType, depth, r3d=r3d2, lfill=lfill, &
                        stat=hit_error)
       IF(test_pass)test_pass = .NOT. hit_error
    END IF
    IF( PRESENT(r3d3) )THEN
       CALL array_check(gridType, testType, depth, r3d=r3d3, lfill=lfill, &
                        stat=hit_error)
       IF(test_pass)test_pass = .NOT. hit_error
    END IF
    IF( PRESENT(r2d1) )THEN
       CALL array_check(gridType, testType, depth, r2d=r2d1, lfill=lfill, &
                        stat=hit_error)
       IF(test_pass)test_pass = .NOT. hit_error
    END IF
    IF( PRESENT(r2d2) )THEN
       CALL array_check(gridType, testType, depth, r2d=r2d2, lfill=lfill, &
                        stat=hit_error)
       IF(test_pass)test_pass = .NOT. hit_error
    END IF
    IF( PRESENT(i2d1) )THEN
       CALL array_check(gridType, testType, depth, i2d=i2d1, lfill=lfill, &
                        stat=hit_error)
       IF(test_pass)test_pass = .NOT. hit_error
    END IF
    IF( PRESENT(i3d1) )THEN
       CALL array_check(gridType, testType, depth, i3d=i3d1, lfill=lfill, &
                        stat=hit_error)
       IF(test_pass)test_pass = .NOT. hit_error
    END IF
    IF( PRESENT(i3d2) )THEN
       CALL array_check(gridType, testType, depth, i3d=i3d2, lfill=lfill, &
                        stat=hit_error)
       IF(test_pass)test_pass = .NOT. hit_error
    END IF

    IF(test_pass)THEN
       WRITE(LOG_UNIT,FMT="('PASS: All arrays OK after exchange.')")
    ELSE
       WRITE(LOG_UNIT,FMT="('FAIL: Errors occurred during this exchange (see above).')")
       ! Set module-wide flag to say that test has failed
       test_failed = .TRUE.
    END IF
    WRITE(LOG_UNIT,*)

  END SUBROUTINE exch_test

  !======================================================================

  SUBROUTINE init_test_arrays(testType, &
                              r2d1, r2d2, r3d1, r3d2, r3d3, i2d1, i3d1, i3d2)
    IMPLICIT none
!FTRANS r3d1 :I :I :z
!FTRANS r3d2 :I :I :z
!FTRANS r3d3 :I :I :z
!FTRANS i3d1 :I :I :z
!FTRANS i3d2 :I :I :z
    INTEGER,            INTENT(in) :: testType ! Type of test to do
    REAL (kind=wp), DIMENSION(:,:),   INTENT(inout), OPTIONAL :: r2d1,r2d2
    REAL (kind=wp), DIMENSION(:,:,:), INTENT(inout), OPTIONAL :: r3d1,r3d2,r3d3
    INTEGER, DIMENSION(:,:,:),        INTENT(inout), OPTIONAL :: i3d1,i3d2
    INTEGER, DIMENSION(:,:),          INTENT(inout), OPTIONAL :: i2d1

    ! Initialise array using global coordinates (mig and mjg arrays are
    ! not used here as they've not been set-up yet).

    IF(testType .EQ. GLOBAL_LOCN_TEST)THEN
       IF( PRESENT(r3d1) )THEN
          CALL set_by_global_coords_3dr(r3d1)
       END IF
       IF( PRESENT(r3d2) )THEN
          CALL set_by_global_coords_3dr(r3d2)
       END IF
       IF( PRESENT(r3d3) )THEN
          CALL set_by_global_coords_3dr(r3d3)
       END IF
       IF( PRESENT(r2d1) )THEN
          CALL set_by_global_coords_2dr(r2d1)
       END IF
       IF( PRESENT(r2d2) )THEN
          CALL set_by_global_coords_2dr(r2d2)
       END IF
       IF( PRESENT(i2d1) )THEN
          CALL set_by_global_coords_2di(i2d1)
       END IF
       IF( PRESENT(i3d1) )THEN
          CALL set_by_global_coords_3di(i3d1)
       END IF
       IF( PRESENT(i3d2) )THEN
          CALL set_by_global_coords_3di(i3d2)
       END IF
    ELSE IF(testType .EQ. PE_ID_TEST)THEN
       IF(PRESENT(r3d1) )THEN
          CALL set_by_pe_id_3dr(r3d1)
       END IF
       IF(PRESENT(r3d2) )THEN
          CALL set_by_pe_id_3dr(r3d2)
       END IF
       IF(PRESENT(r3d3) )THEN
          CALL set_by_pe_id_3dr(r3d3)
       END IF
       IF(PRESENT(r2d1) )THEN
          CALL set_by_pe_id_2dr(r2d1)
       END IF
       IF(PRESENT(r2d2) )THEN
          CALL set_by_pe_id_2dr(r2d2)
       END IF
       IF( PRESENT(i2d1) )THEN
          CALL set_by_pe_id_2di(i2d1)
       END IF
       IF( PRESENT(i3d1) )THEN
          CALL set_by_pe_id_3di(i3d1)
       END IF
       IF( PRESENT(i3d2) )THEN
          CALL set_by_pe_id_3di(i3d2)
       END IF
    END IF

  END SUBROUTINE init_test_arrays

  !======================================================================

  SUBROUTINE set_by_pe_id_3dr(r3d)
    USE par_kind,     ONLY: wp
    USE dom_oce,      ONLY: narea, nlci, nlcj
    USE mapcomm_mod,  ONLY: ilbext, iubext, jlbext, jubext
    IMPLICIT none
!FTRANS r3d :I :I :z
    ! Arguments
    REAL (kind=wp), DIMENSION(:,:,:), INTENT(inout) :: r3d

    ! Initialise array with the PE id
    r3d(:, :, :) = narea

    ! Set halos so that they contain negative values. If we do not
    ! have cyclic boundary conditions then there are no halos on the E and W
    ! boundaries of the model domain.
    IF((.NOT. ilbext) .OR. (ilbext .AND. cyclic_bc) ) r3d(1, :, :)   = -narea
    IF((.NOT. iubext) .OR. (iubext .AND. cyclic_bc) )THEN
       r3d(nlci:,:,:) = -narea
    END IF

    ! No halos at top and bottom of global domain
    IF(.not.jlbext)r3d(:,1,:) = -narea
    IF(.not.jubext)r3d(:,nlcj:,:) = -narea

  END SUBROUTINE set_by_pe_id_3dr

  !======================================================================

  SUBROUTINE set_by_pe_id_2dr(r2d)
    USE par_kind,     ONLY: wp
    USE dom_oce,      ONLY: narea, nlci, nlcj
    USE mapcomm_mod,  ONLY: ilbext, iubext, jlbext, jubext
    IMPLICIT none
    ! Arguments
    REAL (kind=wp), DIMENSION(:,:), INTENT(inout) :: r2d

    ! Initialise array with the PE id
    r2d(:, :) = narea

    ! Set halos so that they contain negative values. If we do not
    ! have cyclic boundary conditions then there are no halos on the E and W
    ! boundaries of the model domain.
    IF((.NOT. ilbext) .OR. (ilbext .AND. cyclic_bc) ) r2d(1, :)    = -narea
    IF((.NOT. iubext) .OR. (iubext .AND. cyclic_bc) )THEN
       r2d(nlci:,:) = -narea
    END IF

    ! No halos at top and bottom of global domain
    IF(.not.jlbext)r2d(:,1)     = -narea
    IF(.not.jubext)r2d(:,nlcj:) = -narea

  END SUBROUTINE set_by_pe_id_2dr

  !======================================================================

  SUBROUTINE set_by_pe_id_3di(i3d)
    USE dom_oce,      ONLY: narea, nlci, nlcj
    USE mapcomm_mod,  ONLY: ilbext, iubext, jlbext, jubext
    IMPLICIT none
!FTRANS i3d :I :I :z
    ! Arguments
    INTEGER, DIMENSION(:,:,:), INTENT(inout) :: i3d

    ! Initialise array with the PE id
    i3d(:, :, :) = narea

    ! Set halos so that they contain negative values. If we do not
    ! have cyclic boundary conditions then there are no halos on the E and W
    ! boundaries of the model domain.
    IF((.NOT. ilbext) .OR. (ilbext .AND. cyclic_bc) ) i3d(1, :, :)   = -narea
    IF((.NOT. iubext) .OR. (iubext .AND. cyclic_bc) )THEN
       i3d(nlci:,:,:) = -narea
    END IF

    ! No halos at top and bottom of global domain
    IF(.not.jlbext)i3d(:,1,:) = -narea
    IF(.not.jubext)i3d(:,nlcj:,:) = -narea

  END SUBROUTINE set_by_pe_id_3di

  !======================================================================

  SUBROUTINE set_by_pe_id_2di(i2d)
    USE dom_oce,      ONLY: narea, nlci, nlcj
    USE mapcomm_mod,  ONLY: ilbext, iubext, jlbext, jubext
    IMPLICIT none
    ! Arguments
    INTEGER, DIMENSION(:,:), INTENT(inout) :: i2d
    !!-------------------------------------------------------------------

    ! Initialise array with the PE id
    i2d(:, :) = narea

    ! Set halos so that they contain negative values. If we do not
    ! have cyclic boundary conditions then there are no halos on the E and W
    ! boundaries of the model domain.
    IF((.NOT. ilbext) .OR. (ilbext .AND. cyclic_bc) )i2d(1, :)   = -narea
    IF((.NOT. iubext) .OR. (iubext .AND. cyclic_bc) )THEN
       i2d(nlci:,:) = -narea
    END IF

    ! No halos at top and bottom of global domain
    IF(.not.jlbext)i2d(:,1) = -narea
    IF(.not.jubext)i2d(:,nlcj:) = -narea

  END SUBROUTINE set_by_pe_id_2di

  !======================================================================

  SUBROUTINE set_by_global_coords_3dr(r3d)
    USE par_kind,     ONLY: wp
    USE par_oce,      ONLY: jpi, jpj, jpk, jpreci
    USE dom_oce,      ONLY: nlci, nlcj, nldi, nlei, &
                            nldj, nlej, narea
    USE mapcomm_mod,  ONLY: ilbext, iubext, jlbext, jubext
    IMPLICIT none
!FTRANS r3d :I :I :z
    ! Arguments
    REAL (kind=wp), DIMENSION(:,:,:), INTENT(inout) :: r3d
    ! Locals
    INTEGER :: ik, ij, ii
    !!-------------------------------------------------------------------

    ! Initialise array using global coordinates (mig and mjg arrays are
    ! not used here as they've not been set-up yet).
    ! A correct exchange process (but without north-fold) won't change
    ! these values
#if defined key_z_first
    DO ij=1,jpj,1
       DO ii=1,jpi,1
          DO ik=1,jpk,1
#else
    DO ik=1,jpk,1
       DO ij=1,jpj,1
          DO ii=1,jpi,1
#endif
             r3d(ii,ij,ik) = REAL( gcoords_to_int(ii,ij,ik) )
          END DO
       END DO
    END DO

    ! Set halos so that they contain negative values. If we do not
    ! have cyclic boundary conditions then there are no halos on the E and W
    ! boundaries of the model domain.
    IF( (.NOT. ilbext) .OR. (ilbext .AND. cyclic_bc) )r3d(1, :, :)   = -1.0_wp*r3d(nldi, :, :)
    IF( (.NOT. iubext) .OR. (iubext .AND. cyclic_bc) )THEN
       DO ii=nlci,jpi,1
          r3d(ii,:,:) = -1.0_wp*r3d(nlei, :, :)
       END DO
    END IF

    ! No halos at top and bottom of global domain
    IF(.NOT.jlbext)r3d(:,1,:) = -1.0_wp*r3d(:,2,:)
    IF(.NOT.jubext)THEN
       DO ij=nlcj, jpj, 1
          r3d(:,ij,:) = -1.0_wp*r3d(:,nlej,:)
       END DO
    END IF

  END SUBROUTINE set_by_global_coords_3dr

  !======================================================================

  SUBROUTINE set_by_global_coords_2dr(r2d)
    USE par_kind,     ONLY: wp
    USE par_oce,      ONLY: jpi, jpj, jpk, jpreci
    USE dom_oce,      ONLY: nlci, nlcj, nldi, nlei, &
                            nldj, nlej
    USE mapcomm_mod,  ONLY: ilbext, iubext, jlbext, jubext
    IMPLICIT none
    REAL (kind=wp), DIMENSION(:,:), INTENT(inout) :: r2d
    ! Locals
    INTEGER :: ij, ii

    ! Initialise array using global coordinates (mig and mjg arrays are
    ! not used here as they've not been set-up yet).
    ! A correct exchange process (but without north-fold) won't change
    ! these values
    DO ij=nldj,jpj,1
       DO ii=nldi,jpi,1
          r2d(ii,ij) = REAL( gcoords_to_int(ii,ij) )
       END DO
    END DO

    ! Set halos so that they contain negative values. If we do not
    ! have cyclic boundary conditions then there are no halos on the E and W
    ! boundaries of the model domain.
    IF((.NOT. ilbext) .OR. (ilbext .AND. cyclic_bc) )r2d(1, :)   = -1.0*r2d(nldi, :)
    IF((.NOT. iubext) .OR. (iubext .AND. cyclic_bc) )THEN
       DO ii=nlci,jpi,1
          r2d(ii,:) = -1.0*r2d(nlei, :)
       END DO
    END IF
    ! No halos at top and bottom of global domain
    IF(.NOT.jlbext)r2d(:,1) = -1.0_wp*r2d(:,2)
    IF(.NOT.jubext)THEN
       DO ij=nlcj, jpj, 1
          r2d(:,ij) = -1.0*r2d(:,nlej)
       END DO
    END IF

  END SUBROUTINE set_by_global_coords_2dr

  !======================================================================

  SUBROUTINE compare_arrays2d(r2d1, r2d2)
    USE par_kind,     ONLY: wp
    IMPLICIT none
    REAL (kind=wp), DIMENSION(:,:),   INTENT(in) :: r2d1,r2d2

    !WRITE(*,*) 'compare_arrays2d: IMPLEMENT ME!'

  END SUBROUTINE compare_arrays2d

  !======================================================================

  SUBROUTINE compare_arrays3d(r3d1, r3d2)
    USE par_kind,     ONLY: wp
    IMPLICIT none
    REAL (kind=wp), DIMENSION(:,:,:), INTENT(in) :: r3d1,r3d2

    !WRITE(*,*) 'compare_arrays3d: IMPLEMENT ME!'

  END SUBROUTINE compare_arrays3d

  !======================================================================

  SUBROUTINE exch_timing(gridType, &
                         r2d1, r2d2, &
                         r3d1, r3d2, r3d3, &
                         i2d1, i2d2, i3d1, i3d2, isgn, lfill)
    USE par_kind,     ONLY: wp
    USE par_oce,      ONLY: jpreci
    USE mapcomm_mod,  ONLY: Iminus, Iplus, Jminus, Jplus, NONE
    USE mpi,          ONLY: MPI_COMM_WORLD, MPI_Wtime
    IMPLICIT none
    ! Arguments
    CHARACTER(LEN=1),   INTENT(in) :: gridType ! Grid on which test being done
    REAL (kind=wp), DIMENSION(:,:),   INTENT(inout), OPTIONAL :: r2d1,r2d2
    REAL (kind=wp), DIMENSION(:,:,:), INTENT(inout), OPTIONAL :: r3d1,r3d2,r3d3
    INTEGER, DIMENSION(:,:,:),        INTENT(inout), OPTIONAL :: i3d1,i3d2
    INTEGER, DIMENSION(:,:),          INTENT(inout), OPTIONAL :: i2d1,i2d2
    INTEGER,                             INTENT(in), OPTIONAL :: isgn
    LOGICAL,                             INTENT(in), OPTIONAL :: lfill
    ! Local vars
    INTEGER :: iloop, narrays
    REAL    :: t1, t2, ttaken, speed
    INTEGER :: nr3darrays, nr2darrays, ni3darrays, ni2darrays
    LOGICAL :: lfilllocal
    CHARACTER(LEN=200) :: fmtString

    IF(PRESENT(lfill))THEN
       lfilllocal = lfill
    ELSE
       lfilllocal = .FALSE.
    END IF

    narrays = 0
    nr3darrays = 0
    nr2darrays = 0
    ni3darrays = 0
    ni2darrays = 0

    IF(PRESENT(r3d1) )THEN
       CALL set_by_pe_id_3dr(r3d1)
       nr3darrays = nr3darrays + 1
    END IF
    IF(PRESENT(r3d2) )THEN
       CALL set_by_pe_id_3dr(r3d2)
       nr3darrays = nr3darrays + 1
    END IF
    IF(PRESENT(r3d3) )THEN
       CALL set_by_pe_id_3dr(r3d3)
       nr3darrays = nr3darrays + 1
    END IF
    IF(PRESENT(r2d1) )THEN
       CALL set_by_pe_id_2dr(r2d1)
       nr2darrays = nr2darrays +1
    END IF
    IF(PRESENT(r2d2) )THEN
       CALL set_by_pe_id_2dr(r2d2)
       nr2darrays = nr2darrays +1
    END IF
    IF( PRESENT(i2d1) )THEN
       CALL set_by_pe_id_2di(i2d1)
       ni2darrays = ni2darrays +1
    END IF
    IF( PRESENT(i2d2) )THEN
       CALL set_by_pe_id_2di(i2d2)
       ni2darrays = ni2darrays +1
    END IF
    IF( PRESENT(i3d1) )THEN
       CALL set_by_pe_id_3di(i3d1)
       ni3darrays = ni3darrays +1
    END IF
    IF( PRESENT(i3d2) )THEN
       CALL set_by_pe_id_3di(i3d2)
       ni3darrays = ni3darrays +1
    END IF

    narrays = nr3darrays + nr2darrays + ni3darrays + ni2darrays

    IF(use_exch_list)THEN

       CALL time_packed_exch(ttaken, gridType,                &
                             r2d1=r2d1, r2d2=r2d2,            &
                             r3d1=r3d1, r3d2=r3d2, r3d3=r3d3, &
                             i2d1=i2d1, i2d2=i2d2,            &
                             i3d1=i3d1, i3d2=i3d2,            &
                             isgn=isgn, lfill=lfill)

#ifdef key_mpp_mpi
       IF (ttaken > 0.0D0) THEN
          speed = dble(loops_to_time*narrays)/ttaken
       ELSE
          speed = 0.0D0
       END IF

       ! This format string exceeds the maximum fortran line length and must
       ! therefore be broken up. However, using the concatenation operator
       ! 'in place' within the write statement has nasty side effects on the
       ! Cray XE with the cray compiler (maybe stack related?). Therefore we 
       ! use fmtString as a temporary variable to construct the string.
       fmtString = ""
       fmtString = "(I6,'*(',I1,' 2Di,',I1,' 3Di,',I1,' 2Dr,',I1,' 3Dr) "// &
                   "packd exchs at ',A,' in ',F10.6,'s, ',F10.1,' arr ex/s, lfill = ',L1)"
       WRITE (LOG_UNIT, TRIM(fmtString)) &
        loops_to_time,ni2darrays,ni3darrays,nr2darrays,nr3darrays,gridType,&
        ttaken,speed, lfilllocal
#endif
    ELSE
       WRITE (LOG_UNIT, "('Packed exchanges not timed because exchange lists switched off in exchtestmod.F90')")
    END IF

    CALL time_indiv_exch(ttaken, gridType,                &
                         r2d1=r2d1, r2d2=r2d2,            &
                         r3d1=r3d1, r3d2=r3d2, r3d3=r3d3, &
                         i2d1=i2d1, i2d2=i2d2,            &
                         i3d1=i3d1, i3d2=i3d2,            &
                         isgn=isgn, lfill=lfill)

#ifdef key_mpp_mpi
    IF (ttaken > 0.0D0) THEN
       speed = dble(loops_to_time*narrays)/ttaken
    ELSE
       speed = 0.0D0
    END IF

    ! See comment on first use of fmtString above.
    fmtString = ""
    fmtString = "(I6,'*(',I1,' 2Di,',I1,' 3Di,',I1,' 2Dr,',I1,' 3Dr) "// &
                "indiv. exchanges at ',A,' in ',F10.6,'s, ',F10.1,' array exch/s, lfill = ',L1)"
    WRITE (LOG_UNIT, TRIM(fmtString)) &
        loops_to_time,ni2darrays,ni3darrays,nr2darrays,nr3darrays,gridType,&
        ttaken,speed, lfilllocal
!!$    WRITE (LOG_UNIT,"(I6,'*(',I1,' 2Di,',I1,' 3Di,',I1,' 2Dr,',I1,' 3Dr) "// &
!!$"indiv. exchanges at ',A,' in ',F10.6,'s, ',F10.1,' array exch/s, lfill = ',L1)") &
!!$        loops_to_time,ni2darrays,ni3darrays,nr2darrays,nr3darrays,gridType,&
!!$        ttaken,speed, lfilllocal
#endif

  END SUBROUTINE exch_timing

  !======================================================================

  SUBROUTINE time_packed_exch(ttaken, gridType, &
                              r2d1, r2d2, &
                              r3d1, r3d2, r3d3, &
                              i2d1, i2d2, i3d1, i3d2, &
                              isgn, lfill)
    USE par_kind,     ONLY: wp
    USE par_oce,      ONLY: jpreci
    USE mapcomm_mod,  ONLY: Iminus, Iplus, Jminus, Jplus, NONE
    USE exchmod,      ONLY: add_exch, bound_exch_list
    USE mpi,          ONLY: MPI_COMM_WORLD, MPI_Wtime
    IMPLICIT none
    ! Arguments
    REAL (kind=wp),    INTENT(out) :: ttaken
    CHARACTER(LEN=1),   INTENT(in) :: gridType ! Grid on which test being done
    REAL (kind=wp), DIMENSION(:,:),   INTENT(inout), OPTIONAL :: r2d1,r2d2
    REAL (kind=wp), DIMENSION(:,:,:), INTENT(inout), OPTIONAL :: r3d1,r3d2,r3d3
    INTEGER, DIMENSION(:,:,:),        INTENT(inout), OPTIONAL :: i3d1,i3d2
    INTEGER, DIMENSION(:,:),          INTENT(inout), OPTIONAL :: i2d1,i2d2
    INTEGER,                             INTENT(in), OPTIONAL :: isgn
    LOGICAL,                             INTENT(in), OPTIONAL :: lfill
    ! Local vars
    REAL    :: t1, t2
    INTEGER :: iloop, ierror

#if ! defined key_mpp_rkpart
    ttaken = 0.0
    RETURN
#endif

#ifdef key_mpp_mpi
    CALL MPI_barrier(MPI_comm_world, ierror)
    t1 = MPI_wtime()
#endif

    DO iloop = 1, loops_to_time, 1

       IF(PRESENT(r3d1))THEN
          CALL add_exch(jpreci, gridType, Iminus,Jminus,Iplus,Jplus, r3d=r3d1, &
                        isgn=isgn, lfill=lfill)
       END IF
       IF(PRESENT(r3d2))THEN
          CALL add_exch(jpreci, gridType, Iminus,Jminus,Iplus,Jplus, r3d=r3d2, &
                        isgn=isgn, lfill=lfill)
       END IF
       IF(PRESENT(r3d3))THEN
          CALL add_exch(jpreci, gridType, Iminus,Jminus,Iplus,Jplus, r3d=r3d3, &
                        isgn=isgn, lfill=lfill)
       END IF
       IF(PRESENT(r2d1))THEN
          CALL add_exch(jpreci, gridType, Iminus,Jminus,Iplus,Jplus, r2d=r2d1, &
                        isgn=isgn, lfill=lfill)
       END IF
       IF(PRESENT(r2d2))THEN
          CALL add_exch(jpreci, gridType, Iminus,Jminus,Iplus,Jplus, r2d=r2d2, &
                        isgn=isgn, lfill=lfill)
       END IF
       IF( PRESENT(i2d1) )THEN
          CALL add_exch(jpreci, gridType, Iminus,Jminus,Iplus,Jplus, i2d=i2d1, &
                        isgn=isgn, lfill=lfill)
       END IF
       IF( PRESENT(i2d2) )THEN
          CALL add_exch(jpreci, gridType, Iminus,Jminus,Iplus,Jplus, i2d=i2d2, &
                        isgn=isgn, lfill=lfill)
       END IF
       IF( PRESENT(i3d1) )THEN
          CALL add_exch(jpreci, gridType, Iminus,Jminus,Iplus,Jplus, i3d=i3d1, &
                        isgn=isgn, lfill=lfill)
       END IF
       IF( PRESENT(i3d2) )THEN
          CALL add_exch(jpreci, gridType, Iminus,Jminus,Iplus,Jplus, i3d=i3d2, &
                        isgn=isgn, lfill=lfill)
       END IF
       
       CALL bound_exch_list()

    END DO

#ifdef key_mpp_mpi
    t2 = MPI_wtime()
    ttaken = t2-t1
#endif

  END SUBROUTINE time_packed_exch

  !======================================================================

  SUBROUTINE time_indiv_exch(ttaken, gridType, &
                             r2d1, r2d2, &
                             r3d1, r3d2, r3d3, &
                             i2d1, i2d2, i3d1, i3d2, &
                             isgn, lfill)
    USE par_kind,     ONLY: wp
    USE par_oce,      ONLY: jpreci
    USE mapcomm_mod,  ONLY: Iminus, Iplus, Jminus, Jplus, NONE
    USE exchmod,      ONLY: bound_exch
    USE lbclnk,       ONLY: lbc_lnk
    USE mpi,          ONLY: MPI_COMM_WORLD, MPI_Wtime
    IMPLICIT none
    ! Arguments
    REAL (kind=wp),    INTENT(out) :: ttaken
    CHARACTER(LEN=1),   INTENT(in) :: gridType ! Grid on which test being done
    REAL (kind=wp), DIMENSION(:,:),   INTENT(inout), OPTIONAL :: r2d1,r2d2
    REAL (kind=wp), DIMENSION(:,:,:), INTENT(inout), OPTIONAL :: r3d1,r3d2,r3d3
    INTEGER, DIMENSION(:,:,:),        INTENT(inout), OPTIONAL :: i3d1,i3d2
    INTEGER, DIMENSION(:,:),          INTENT(inout), OPTIONAL :: i2d1,i2d2
    INTEGER,                             INTENT(in), OPTIONAL :: isgn
    LOGICAL,                             INTENT(in), OPTIONAL :: lfill
    ! Local vars
    REAL    :: t1, t2, psgn
    INTEGER :: iloop, ierror

#ifdef key_mpp_mpi
    CALL MPI_barrier(MPI_comm_world, ierror)
    t1 = MPI_wtime()
#endif

    indiv_timing: DO iloop = 1, loops_to_time, 1

       IF(PRESENT(r3d1))THEN
#if defined key_mpp_rkpart
          CALL bound_exch(r3d1, jpreci, jpreci,         &
                          Iplus, Iminus, Jplus, Jminus, &
                          gridType, isgn=isgn, lfill=lfill)
#else
          ! psgn is NOT an optional argument to lbc_lnk
          psgn = 1.0_wp
          IF( PRESENT(isgn) )psgn=REAL(isgn)
          CALL lbc_lnk(r3d1, gridType, psgn)
#endif
       END IF
       IF(PRESENT(r3d2))THEN
#if defined key_mpp_rkpart
          CALL bound_exch(r3d2, jpreci, jpreci,         &
                          Iplus, Iminus, Jplus, Jminus, &
                          gridType, isgn=isgn, lfill=lfill)
#endif
       END IF
       IF(PRESENT(r3d3))THEN
#if defined key_mpp_rkpart
          CALL bound_exch(r3d3, jpreci, jpreci,         &
                          Iplus, Iminus, Jplus, Jminus, &
                          gridType, isgn=isgn, lfill=lfill)
#endif
       END IF
       IF(PRESENT(r2d1))THEN
#if defined key_mpp_rkpart
          CALL bound_exch(r2d1, jpreci, jpreci,         &
                          Iplus, Iminus, Jplus, Jminus, &
                          gridType, isgn=isgn, lfill=lfill)
#else
          ! psgn is NOT an optional argument to lbc_lnk
          psgn = 1.0_wp
          IF( PRESENT(isgn) )psgn=REAL(isgn)
          CALL lbc_lnk(r2d1, gridType, psgn)
#endif
       END IF
       IF(PRESENT(r2d2))THEN
#if defined key_mpp_rkpart
          CALL bound_exch(r2d2, jpreci, jpreci,         &
                          Iplus, Iminus, Jplus, Jminus, &
                          gridType, isgn=isgn, lfill=lfill)
#else
          ! psgn is NOT an optional argument to lbc_lnk
          psgn = 1.0_wp
          IF( PRESENT(isgn) )psgn=REAL(isgn)
          CALL lbc_lnk(r2d2, gridType, psgn)
#endif
       END IF
       IF( PRESENT(i2d1) )THEN
#if defined key_mpp_rkpart
          CALL bound_exch(i2d1, jpreci, jpreci,         &
                          Iplus, Iminus, Jplus, Jminus, &
                          gridType, isgn=isgn, lfill=lfill)
#endif
       END IF
       IF( PRESENT(i2d2) )THEN
#if defined key_mpp_rkpart
          CALL bound_exch(i2d2, jpreci, jpreci,         &
                          Iplus, Iminus, Jplus, Jminus, &
                          gridType, isgn=isgn, lfill=lfill)
#endif
       END IF
       IF( PRESENT(i3d1) )THEN
#if defined key_mpp_rkpart
          CALL bound_exch(i3d1, jpreci, jpreci,         &
                          Iplus, Iminus, Jplus, Jminus, &
                          gridType, isgn=isgn, lfill=lfill)
#endif
       END IF
       IF( PRESENT(i3d2) )THEN
#if defined key_mpp_rkpart
          CALL bound_exch(i3d2, jpreci, jpreci,         &
                          Iplus, Iminus, Jplus, Jminus, &
                          gridType, isgn=isgn, lfill=lfill)
#endif
       END IF
       
    END DO indiv_timing

#ifdef key_mpp_mpi
    t2 = MPI_wtime()
    ttaken = t2-t1
#endif

  END SUBROUTINE time_indiv_exch

  !======================================================================

  SUBROUTINE set_by_global_coords_3di(i3d)
    USE par_kind,     ONLY: wp
    USE par_oce,      ONLY: jpi, jpj, jpk, jpreci
    USE dom_oce,      ONLY: nlci, nlcj, nldi, nlei, &
                            nldj, nlej
    USE mapcomm_mod,  ONLY: ilbext, iubext, jlbext, jubext
    IMPLICIT none
!FTRANS i3d :I :I :z
    ! Arguments
    INTEGER, DIMENSION(:,:,:), INTENT(inout) :: i3d
    ! Locals
    INTEGER :: ik, ij, ii

    ! Initialise array using global coordinates (mig and mjg arrays are
    ! not used here as they've not been set-up yet).
    ! A correct exchange process (but without north-fold) won't change
    ! these values
#if defined key_z_first
    DO ij=nldj,jpj,1
       DO ii=nldi,jpi,1
          DO ik=1,jpk,1
#else
    DO ik=1,jpk,1
       DO ij=nldj,jpj,1
          DO ii=nldi,jpi,1
#endif
             i3d(ii,ij,ik) = gcoords_to_int(ii,ij,ik)
          END DO
       END DO
    END DO

    ! Set halos so that they contain negative values. If we do not
    ! have cyclic boundary conditions then there are no halos on the E and W
    ! boundaries of the model domain.
    IF((.NOT. ilbext) .OR. (ilbext .AND. cyclic_bc) )i3d(1, :, :)   = -1*i3d(nldi, :, :)
    IF((.NOT. iubext) .OR. (iubext .AND. cyclic_bc) )THEN
       DO ii=nlci,jpi,1
          i3d(ii,:,:) = -1*i3d(nlei, :, :)
       END DO
    END IF

    ! No halos at top and bottom of global domain
    IF(.not.jlbext)i3d(:,1,:) = -1*i3d(:,2,:)
    IF(.not.jubext)THEN
       DO ij=nlcj, jpj, 1
          i3d(:,ij,:) = -1*i3d(:,nlej,:)
       END DO
    END IF

  END SUBROUTINE set_by_global_coords_3di

  !======================================================================

  SUBROUTINE set_by_global_coords_2di(i2d)
    USE par_kind,     ONLY: wp
    USE par_oce,      ONLY: jpi, jpj, jpk, jpreci
    USE dom_oce,      ONLY: nlci, nlcj, nldi, nlei, &
                            nldj, nlej
    USE mapcomm_mod,  ONLY: ilbext, iubext, jlbext, jubext
    IMPLICIT none
    ! Arguments
    INTEGER, DIMENSION(:,:), INTENT(inout) :: i2d
    ! Locals
    INTEGER :: ij, ii

    ! Initialise array using global coordinates (mig and mjg arrays are
    ! not used here as they've not been set-up yet).
    ! A correct exchange process (but without north-fold) won't change
    ! these values
    DO ij=nldj,jpj,1
       DO ii=nldi,jpi,1
          i2d(ii,ij) = gcoords_to_int(ii,ij)
       END DO
    END DO

    ! Set halos so that they contain negative values
    i2d(1, :)   = -1*i2d(nldi, :)
    DO ii=nlci,jpi,1
       i2d(ii,:) = -1*i2d(nlei, :)
    END DO
    ! No halos at top and bottom of global domain
    IF(.not.jlbext)i2d(:,1) = -1*i2d(:,2)
    IF(.not.jubext)THEN
       DO ij=nlcj, jpj, 1
          i2d(:,ij) = -1*i2d(:,nlej)
       END DO
    END IF

  END SUBROUTINE set_by_global_coords_2di

  !======================================================================

  SUBROUTINE array_check(gridType, testType, depth, r2d, r3d, i2d, i3d, &
                         lfill, stat)
    USE par_kind,     ONLY: wp
    USE par_oce,      ONLY: jpi, jpj, jpk, jpreci, jpiglo
    USE dom_oce,      ONLY: nlci, nldi, nlei, nldj, nlej, nimpp, njmpp, narea
    USE mapcomm_mod,  ONLY: jlbext, jubext, ilbext, iubext
    USE exchmod,      ONLY: num_nfold_rows
    IMPLICIT none
!FTRANS r3d :I :I :z
!FTRANS i3d :I :I :z
    ! Arguments
    CHARACTER(LEN=1), INTENT(in) :: gridType
    INTEGER,          INTENT(in) :: testType
    INTEGER,        DIMENSION(:,:),   INTENT(in)           :: depth ! Global land mask
    REAL (kind=wp), DIMENSION(:,:),   INTENT(in), OPTIONAL :: r2d
    REAL (kind=wp), DIMENSION(:,:,:), INTENT(in), OPTIONAL :: r3d
    INTEGER,        DIMENSION(:,:),   INTENT(in), OPTIONAL :: i2d
    INTEGER,        DIMENSION(:,:,:), INTENT(in), OPTIONAL :: i3d
    LOGICAL,                          INTENT(in), OPTIONAL :: lfill
    LOGICAL,                          INTENT(out)          :: stat
    ! Locals
    INTEGER :: ik, ij, ii, ipt
    LOGICAL :: hit_error, local_lfill
    INTEGER :: gVal, jstart, jstop, istart, istop
    REAL (kind=wp), PARAMETER :: TOL_ZERO = 1.0e-6

    hit_error = .FALSE.

    IF(PRESENT(lfill))THEN
       local_lfill = lfill
    ELSE
       local_lfill = .FALSE.
    END IF

    WRITE (LOG_UNIT, &
           FMT="('array_check: grid = ',A,' test = ',I1, ' lfill = ',L1)") &
           gridType, testType, local_lfill

    ! Exclude halos either side of GLOBAL domain because of cyclic
    ! B.C.'s
    istart = 1
    ! Halos at either side of simulation domain are zeroed if the domain
    ! is closed and lfill is false.
    IF( ilbext .AND. (.NOT. local_lfill) .AND. (.NOT. cyclic_bc) )istart = 2
    !IF(ilbext .AND. cyclic_bc)istart = 2
    istop = nlci
    IF( iubext .AND. (.NOT. local_lfill) .AND. (.NOT. cyclic_bc) )istop = nlci-1
!    IF(iubext .AND. cyclic_bc)istop = nlci - 1

    ! No halos at top and bottom of GLOBAL domain
    jstart = nldj
    IF(jlbext)THEN
       IF(local_lfill)THEN
          jstart = 1
       ELSE
          jstart = 2 ! Bottom row of global domain is set to zero
                     ! when lfill is false.
       END IF
    END IF

    jstop = nlej
    IF(jubext)THEN
       IF(local_lfill)THEN
          jstop = nlej
       ELSE
          ! If lfill isn't true then the north-fold condition is applied
          ! and our simple test will fail so exclude any affected rows
          jstop = nlej - num_nfold_rows
       END IF
    END IF
       
    ! Now the tests themselves...

    IF( PRESENT(r3d) )THEN

       ! Halo on W edge of model domain is zeroed when lfill is false and
       ! domain is closed (no cyclic b.c.)
       IF( (.NOT. local_lfill) .AND. ilbext .AND. (.NOT. trimmed(widx,narea)) &
            .AND. (.NOT. cyclic_bc) )THEN

#if defined key_z_first
          DO ij = 1, jpj, 1
             DO ii=1,nldi-1,1
                DO ik=1,jpk,1
#else
          DO ik=1,jpk,1
             DO ij = 1, jpj, 1
                DO ii=1,nldi-1,1
#endif
                IF( ABS(depth(ii+nimpp-1,ij+njmpp-1)*r3d(ii,ij,ik)) > TOL_ZERO)THEN
                   hit_error = .TRUE.
                   WRITE (LOG_UNIT,FMT="(I4,': ERROR: ',A1,' pt, r3d(',2(I3,','),I3,') = ',F10.0, ' NOT zero')") &
#if defined key_z_first
                              narea-1, gridType, ik, ii, ij, r3d(ii, 1, ik)
#else
                              narea-1, gridType, ii, ij, ik, r3d(ii, 1, ik)
#endif
                END IF
             END DO
          END DO
          END DO
       END IF

       ! Halo on E edge of model domain is zeroed when lfill is false and
       ! domain is closed (no cyclic b.c.)
       IF( (.NOT. local_lfill) .AND. iubext .AND. (.NOT. trimmed(eidx,narea)) & 
            .AND. (.NOT. cyclic_bc) )THEN

#if defined key_z_first
          DO ij = 1, jpj, 1
             DO ii=nlei+1,jpi,1
                DO ik=1,jpk,1
#else
          DO ik=1,jpk,1
             DO ij = 1, jpj, 1
                DO ii=nlei+1,jpi,1
#endif
                IF( ABS(depth(ii+nimpp-1,ij+njmpp-1)*r3d(ii,ij,ik)) > TOL_ZERO )THEN
                   hit_error = .TRUE.
                   WRITE (LOG_UNIT,FMT="(I4,': ERROR: ',A1,' pt, r3d(',2(I3,','),I3,') = ',F10.0, ' NOT zero')") &
#if defined key_z_first
                              narea-1, gridType, ik, ii, ij, r3d(ii, ij, ik)
#else
                              narea-1, gridType, ii, ij, ik, r3d(ii, ij, ik)
#endif
                END IF
             END DO
          END DO
          END DO
       END IF

       ! Our code initialises halo regions to negative values. After a 
       ! successful halo swap there should be no negative values anywhere.
#if defined key_z_first
       DO ij = jstart, jstop, 1
          DO ii=1,jpi,1
             DO ik = 1, jpk, 1
#else
       DO ik = 1, jpk, 1
          DO ij = jstart, jstop, 1
             DO ii=1,jpi,1
#endif
                IF(ik <= pmaxdepth(ii+nimpp-1,ij+njmpp-1))THEN
                   IF( depth(ii+nimpp-1,ij+njmpp-1)*r3d(ii,ij,ik) < 0.0 )THEN
                      hit_error = .TRUE.
                      WRITE(LOG_UNIT, &
                            FMT="(I4,': ERROR: ',A1,' pt, r3d(',2(I3,','),I3,') = ',F10.0)") &
#if defined key_z_first
                           narea-1, gridType, ik, ii, ij, r3d(ii, ij, ik)
#else
                           narea-1, gridType, ii, ij, ik, r3d(ii, ij, ik)
#endif
                   END IF
                END IF
             END DO
          END DO
       END DO

       IF(testType .eq. GLOBAL_LOCN_TEST)THEN

          WRITE(LOG_UNIT, *) 'istart, istop = ',istart, istop
          WRITE(LOG_UNIT, *) 'jstart, jstop = ',jstart, jstop
 
          DO ik=1,jpk,1
             DO ij=jstart,jstop,1
                DO ii=istart,istop,1

                   IF(ik <= pmaxdepth(ii+nimpp-1,ij+njmpp-1))THEN

                      ipt = ii+nimpp-1
                      ! Treat halo regions on E/W edges of global domain
                      ! with care when cyclic boundary conditions are
                      ! enabled.
                      IF(cyclic_bc)THEN
                         IF( (ii+nimpp-1) == jpiglo )THEN
                            ! Eastern edge of global domain - this halo
                            ! should therefore contain values from the
                            ! first non-halo column on the Western edge.
                            ipt = 2
                         ELSE IF( (ii+nimpp-1) == 1 )THEN
                            ! Western edge of global domain - this halo
                            ! should therefore contain values from the 
                            ! last non-halo column on the Eastern edge.
                            ipt = jpiglo-1
                         END IF
                      END IF

                      gval = gcoords_to_int(ipt, (ij+njmpp-1), ik, &
                                            are_global=.TRUE.)

                      ! depth is the mask for the whole simulation domain so
                      ! must convert from local to domain coordinates
                      IF( (depth(ii+nimpp-1,ij+njmpp-1) > 0) .AND. &
                          (INT(r3d(ii,ij,ik)) /= gVal) )THEN

                         WRITE(LOG_UNIT, &
                               FMT="(I4,': ERROR: ',A1,' pt, r3d(',2(I3,','),I3,') = ',I10.9,' /= gVal = ',I10.9)") &
#if defined key_z_first
                           narea-1, gridType, ik, ii, ij, &
#else
                           narea-1, gridType, ii, ij, ik, &
#endif
                           INT(r3d(ii, ij, ik)), gVal
                         WRITE (LOG_UNIT,"(I4,': depth(',I3,',',I3,') = ',I2,' bot. level = ',I3)") &
                            narea-1, ii+nimpp-1,ij+njmpp-1, &
                            depth(ii+nimpp-1,ij+njmpp-1),   &
                            pmaxdepth(ii+nimpp-1,ij+njmpp-1)

                         hit_error = .TRUE.
                      END IF
                   END IF

                END DO
             END DO

          END DO

       END IF ! testType == GLOBAL_LOCN_TEST

    END IF ! PRESENT(r3d)

    IF( PRESENT(r2d) )THEN

       DO ij = jstart, jstop, 1
          DO ii=1,jpi,1
             IF( depth(ii,ij)*r2d(ii,ij) < 0.0 )THEN
                hit_error = .TRUE.
                WRITE (LOG_UNIT,FMT="(I4,': ERROR: ',A1,' pt, r2d(',I3,',',I3,') = ',F10.0)") &
                          narea-1, gridType, ii, ij, r2d(ii, ij)
             END IF
          END DO
       END DO

       IF(testType .eq. GLOBAL_LOCN_TEST)THEN

          DO ij=jstart, jstop, 1
             DO ii=istart, istop, 1

                      ipt = ii+nimpp-1
                      ! Treat halo regions on E/W edges of global domain
                      ! with care when cyclic boundary conditions are
                      ! enabled.
                      IF(cyclic_bc)THEN
                         IF( (ii+nimpp-1) == jpiglo )THEN
                            ! Eastern edge of global domain - this halo
                            ! should therefore contain values from the
                            ! first non-halo column on the Western edge.
                            ipt = 2
                         ELSE IF( (ii+nimpp-1) == 1 )THEN
                            ! Western edge of global domain - this halo
                            ! should therefore contain values from the 
                            ! last non-halo column on the Eastern edge.
                            ipt = jpiglo-1
                         END IF
                      END IF

                      gval = gcoords_to_int(ipt, (ij+njmpp-1), &
                                            are_global=.TRUE.)

!                gval = gcoords_to_int(ii,ij)

                IF( (depth(ii+nimpp-1,ij+njmpp-1) > 0) .AND. &
                    (INT(r2d(ii,ij)) /=  gval) )THEN
                   WRITE (LOG_UNIT,FMT="(I4,': ERROR: ',A1,' pt, r2d(',I3,',',I3,') = ',I7.6,' /= gVal = ',I7.6)") &
                        narea-1, gridType, ii, ij, INT(r2d(ii, ij)), gVal
                        hit_error = .TRUE.
                END IF
             END DO
          END DO

       END IF ! testType == GLOBAL_LOCN_TEST

    END IF ! PRESENT(r2d)

    IF( PRESENT(i3d) )THEN

       IF( ANY( MASK=(i3d(1:jpi,jstart:jstop,1:jpk) < 0.0) ) )THEN
          hit_error = .TRUE.
          
          DO ik=1,jpk,1
             DO ij=jstart,jstop,1
                DO ii=1,jpi,1
                   IF( (i3d(ii,ij,ik) < 0.0)  .AND.              &
                       (ik <= pmaxdepth(ii+nimpp-1,ij+njmpp-1)) )THEN

                      WRITE (LOG_UNIT,FMT="(I4,': ERROR: ',A1,' pt, i3d(',2(I3,','),I3,') = ',I10.9)") &
#if defined key_z_first
                              narea-1, gridType, ik, ii, ij, i3d(ii, ij, ik)
#else
                              narea-1, gridType, ii, ij, ik, i3d(ii, ij, ik)
#endif
                   END IF
                END DO
             END DO
          END DO
       END IF ! any array values negative

       IF(testType == GLOBAL_LOCN_TEST)THEN

          WRITE(LOG_UNIT, *) 'istart, istop = ',istart, istop
          WRITE(LOG_UNIT, *) 'jstart, jstop = ',jstart, jstop
 
          DO ik=1,jpk,1

             DO ij=jstart,jstop,1
                DO ii=istart,istop,1
                   gval = gcoords_to_int(ii,ij,ik)
                   IF( (i3d(ii,ij,ik) /= gval) .AND. &
                       (ik <= pmaxdepth(ii+nimpp-1,ij+njmpp-1)) )THEN
                      WRITE (LOG_UNIT,FMT="(I4,': ERROR: ',A1,' pt, i3d(',2(I3,','),I3,') = ',I10.9,' /= gVal = ',I10.9)") &
#if defined key_z_first
                           narea-1, gridType, ik, ii, ij, &
#else
                           narea-1, gridType, ii, ij, ik, &
#endif
                           i3d(ii, ij, ik), gVal
                      hit_error = .TRUE.
                   END IF
                END DO
             END DO

          END DO
       END IF ! testType == GLOBAL_LOCN_TEST

    END IF ! PRESENT(i3d)

    IF( PRESENT(i2d) )THEN

       IF( ANY( MASK=(i2d(1:jpi,jstart:jstop) < 0.0) ) )THEN
          hit_error = .TRUE.
          
          DO ij=jstart,jstop,1
             DO ii=1,jpi,1
                IF(i2d(ii,ij) < 0.0)THEN
                   WRITE (LOG_UNIT,FMT="(I4,': ERROR: ',A1,' pt, i2d(',I3,',',I3,') = ',I10.9)") &
                              narea-1, gridType, ii, ij, i2d(ii, ij)
                END IF
             END DO
          END DO
       END IF ! any array values negative

       IF(testType == GLOBAL_LOCN_TEST)THEN

          WRITE(LOG_UNIT, *) 'istart, istop = ',istart, istop
          WRITE(LOG_UNIT, *) 'jstart, jstop = ',jstart, jstop
 
          DO ij=jstart,jstop,1
             DO ii=istart,istop,1

                gval = gcoords_to_int(ii,ij,ik )
                IF(i2d(ii,ij) /= gval)THEN
                   WRITE (LOG_UNIT,FMT="(I4,': ERROR: ',A1,' pt, i2d(',I3,',',I3,') = ',I10.9,' /= gVal = ',I10.9)") &
                           narea-1, gridType, ii, ij, &
                           i2d(ii, ij), gVal
                   hit_error = .TRUE.
                END IF
             END DO
          END DO

       END IF ! testType == GLOBAL_LOCN_TEST

    END IF ! PRESENT(i2d)

    stat=hit_error

  END SUBROUTINE array_check

  FUNCTION gcoords_to_int(ii, ij, ik, are_global) RESULT(value)
     USE dom_oce,      ONLY: nimpp, njmpp, nldi, nldj
     IMPLICIT None
     ! Convert the specified coordinates in the local domain into global 
     ! coordinates and encode into a single integer number.
     INTEGER, INTENT(in)           :: ii, ij
     INTEGER, INTENT(in), OPTIONAL :: ik
     LOGICAL, INTENT(in), OPTIONAL :: are_global ! Whether input coordinates
                                                 ! are already global rather
                                                 ! than just relative to local
                                                 ! domain
     ! Locals
     INTEGER :: value
     LOGICAL :: lglobal
     !!====================================================================

     lglobal = .FALSE.
     IF( PRESENT(are_global) )lglobal = are_global

     IF(lglobal)THEN
        ! ii and ij are already global coordinates
        value = ii*1000000 + &
                ij*1000
     ELSE
        value = (ii + nimpp - 1)*1000000 + &
                (ij + njmpp - 1)*1000
     END IF

     IF(PRESENT(ik))THEN
        value = value + ik
     END IF

  END FUNCTION gcoords_to_int


END MODULE exchtestmod