Changeset 3837

Timestamp:

2013-03-12T15:55:32+01:00 (11 years ago)

Author:

trackstand2

Message:

Merge of finiss

Location:

branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM

Files:

• ARCH/arch-XE6_HECTOR_CRAY.fcm (added)
• ARCH/arch-ifort_mpi_linux.fcm (added)
• NEMO/LIM_SRC_2/ice_2.F90 (modified) (1 diff)
• NEMO/LIM_SRC_2/limdyn_2.F90 (modified) (1 diff)
• NEMO/LIM_SRC_3/limrhg.F90 (modified) (7 diffs)
• NEMO/OPA_SRC/DOM/closea.F90 (modified) (1 diff)
• NEMO/OPA_SRC/DOM/daymod.F90 (modified) (1 diff)
• NEMO/OPA_SRC/DOM/dom_oce.F90 (modified) (2 diffs)
• NEMO/OPA_SRC/DOM/domain.F90 (modified) (1 diff)
• NEMO/OPA_SRC/DOM/domcfg.F90 (modified) (2 diffs)
• NEMO/OPA_SRC/DOM/domzgr.F90 (modified) (24 diffs)
• NEMO/OPA_SRC/DOM/istate.F90 (modified) (1 diff)
• NEMO/OPA_SRC/DOM/phycst.F90 (modified) (1 diff)
• NEMO/OPA_SRC/DYN/dynnxt.F90 (modified) (2 diffs)
• NEMO/OPA_SRC/DYN/dynspg.F90 (modified) (1 diff)
• NEMO/OPA_SRC/DYN/dynspg_flt.F90 (modified) (2 diffs)
• NEMO/OPA_SRC/DYN/dynzad.F90 (modified) (1 diff)
• NEMO/OPA_SRC/DYN/dynzdf_imp.F90 (modified) (2 diffs)
• NEMO/OPA_SRC/LBC/exchmod.F90 (modified) (121 diffs)
• NEMO/OPA_SRC/LBC/exchtestmod.F90 (modified) (24 diffs)
• NEMO/OPA_SRC/LBC/lbclnk.F90 (modified) (7 diffs)
• NEMO/OPA_SRC/LBC/lib_mpp.F90 (modified) (3 diffs)
• NEMO/OPA_SRC/LBC/mapcomm_mod.F90 (modified) (50 diffs)
• NEMO/OPA_SRC/LBC/mppini.F90 (modified) (2 diffs)
• NEMO/OPA_SRC/LBC/partition_mod.F90 (modified) (35 diffs)
• NEMO/OPA_SRC/LDF/ldfdyn_c3d.h90 (modified) (1 diff)
• NEMO/OPA_SRC/SBC/fldread.F90 (modified) (1 diff)
• NEMO/OPA_SRC/SBC/sbcblk_core.F90 (modified) (3 diffs)
• NEMO/OPA_SRC/SOL/solpcg.F90 (modified) (3 diffs)
• NEMO/OPA_SRC/TRA/traldf_iso.F90 (modified) (2 diffs)
• NEMO/OPA_SRC/ZDF/zdfddm.F90 (modified) (1 diff)
• NEMO/OPA_SRC/arpdebugging.f90 (modified) (10 diffs)
• NEMO/OPA_SRC/nemogcm.F90 (modified) (11 diffs)
• NEMO/OPA_SRC/step.F90 (modified) (3 diffs)
• NEMO/OPA_SRC/timing.F90 (modified) (3 diffs)

branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/LIM_SRC_2/ice_2.F90

@@ -181 +181 @@
       ice_alloc_2 = MAXVAL( ierr )
       !
-      IF( ice_alloc_2 /= 0 )   CALL ctl_warn('ice_alloc_2: failed to allocate arrays')
+      IF( ice_alloc_2 /= 0 )THEN
+         CALL ctl_warn('ice_alloc_2: failed to allocate arrays')
+      ELSE
+         ! Initialise stress tensor to zero
+         stress12_i(:,:) = 0.0_wp
+         stress1_i(:,:)  = 0.0_wp
+         stress2_i(:,:) = 0.0_wp
+      END IF
       !
    END FUNCTION ice_alloc_2

branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/LIM_SRC_2/limdyn_2.F90

@@ -193 +193 @@
             END DO
          END DO
+
          !
       ELSE      ! no ice dynamics : transmit directly the atmospheric stress to the ocean

branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90

@@ -34 +34 @@
    USE dom_ice_2        ! LIM2: ice domain
 #endif
+!   USE arpdebugging, ONLY: dump_array
 
    IMPLICIT NONE
@@ -169 +170 @@
       REAL(wp) ::   zindb         ! ice (1) or not (0)      
       REAL(wp) ::   zdummy        ! dummy argument
+      INTEGER, SAVE :: count = 0 ! For dumping data to disk ARPDBG
       !!-------------------------------------------------------------------
 #if  defined key_lim2 && ! defined key_lim2_vp
@@ -180 +182 @@
      at_i(:,:) = 1. - frld(:,:)
 #endif
+
+      count = count + 1 !  ARPDBG - for dump_array
+
       !
       !------------------------------------------------------------------------------!
@@ -249 +254 @@
       !  v_oce2: ocean v component on v points                        
 
+      !CALL dump_array(count,'tms',tms,withHalos=.TRUE.)
+!      CALL dump_array(count,'e1t',e1t,withHalos=.TRUE.)
+!      CALL dump_array(count,'e2t',e2t,withHalos=.TRUE.)
+      !CALL dump_array(count,'zc1',zc1,withHalos=.TRUE.)
+
       DO jj = k_j1+1, k_jpj-1
          DO ji = fs_2, fs_jpim1
@@ -316 +326 @@
       zs2 (:,:) = stress2_i (:,:)
       zs12(:,:) = stress12_i(:,:)
+
+!      CALL dump_array(count,'u_ice_pre_iter',u_ice,withHalos=.TRUE.)
+!      CALL dump_array(count,'e2u_pre_iter' ,e2u ,withHalos=.TRUE.)
 
       !                                               !----------------------!
@@ -475 +488 @@
             END DO
          END DO
+
          !
          ! Computation of ice velocity
@@ -598 +612 @@
          ENDIF
 
-         !                                                   ! ==================== !
+         !                                                ! ==================== !
       END DO                                              !  end loop over jter  !
       !                                                   ! ==================== !
 
+!      CALL dump_array(count,'u_ice_pre4',u_ice,withHalos=.TRUE.)
       !
       !------------------------------------------------------------------------------!

branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DOM/closea.F90

@@ -32 +32 @@
    PUBLIC clo_bat      ! routine called in domzgr module
 
-   INTEGER, PUBLIC, PARAMETER          ::   jpncs   = 4      !: number of closed sea
-   INTEGER, PUBLIC, DIMENSION(jpncs)   ::   ncstt            !: Type of closed sea
-   INTEGER, PUBLIC, DIMENSION(jpncs)   ::   ncsi1, ncsj1     !: south-west closed sea limits (i,j)
-   INTEGER, PUBLIC, DIMENSION(jpncs)   ::   ncsi2, ncsj2     !: north-east closed sea limits (i,j)
-   INTEGER, PUBLIC, DIMENSION(jpncs)   ::   ncsnr            !: number of point where run-off pours
-   INTEGER, PUBLIC, DIMENSION(jpncs,4) ::   ncsir, ncsjr     !: Location of runoff
-
-   REAL(wp), DIMENSION (jpncs+1)       ::   surf             ! closed sea surface
+   INTEGER, PUBLIC, PARAMETER                ::   jpncs   = 4      !: number of closed sea
+   INTEGER, PUBLIC, SAVE, DIMENSION(jpncs)   ::   ncstt            !: Type of closed sea
+   INTEGER, PUBLIC, SAVE, DIMENSION(jpncs)   ::   ncsi1, ncsj1     !: south-west closed sea limits (i,j)
+   INTEGER, PUBLIC, SAVE, DIMENSION(jpncs)   ::   ncsi2, ncsj2     !: north-east closed sea limits (i,j)
+   INTEGER, PUBLIC, SAVE, DIMENSION(jpncs)   ::   ncsnr            !: number of point where run-off pours
+   INTEGER, PUBLIC, SAVE, DIMENSION(jpncs,4) ::   ncsir, ncsjr     !: Location of runoff
+
+   REAL(wp), SAVE, DIMENSION (jpncs+1)       ::   surf             ! closed sea surface
 
    !! * Control permutation of array indices

branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DOM/daymod.F90

@@ -124 +124 @@
 
       ! control print
-      IF(lwp) WRITE(numout,'(a,i6,a,i2,a,i2,a,i6)')' ==============>> 1/2 time step before the start of the run DATE Y/M/D = ',   &
-           &                   nyear, '/', nmonth, '/', nday, '  nsec_day:', nsec_day, '  nsec_week:', nsec_week
+      IF(lwp) WRITE(numout,"(' ==============>> 1/2 time step before the start of the run DATE Y/M/D = ',i6,'/',i2,'/',i2,'  nsec_day:',i6,'  nsec_week:',I)") nyear, nmonth, nday, nsec_day, nsec_week
 
       ! Up to now, calendar parameters are related to the end of previous run (nit000-1)

branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90

@@ -200 +200 @@
    INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   &
              & mbku, mbkv   !: vertical index of the bottom last U- and W- ocean level
+   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mbkmax       !: Max index of last ocean level on any grid
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bathy        !: ocean depth (meters)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tmask_i      !: interior domain T-point mask
@@ -323 +324 @@
 
       ALLOCATE( mbathy(jpi,jpj) , bathy(jpi,jpj) ,                     &
-         &     tmask_i(jpi,jpj) , bmask(jpi,jpj) ,                     &
+         &     tmask_i(jpi,jpj) , bmask(jpi,jpj) , mbkmax(jpi,jpj),    &
          &     mbkt   (jpi,jpj) , mbku (jpi,jpj) , mbkv(jpi,jpj) , STAT=ierr(9) )
 

branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90

@@ -40 +40 @@
 
    PUBLIC   dom_init   ! called by opa.F90
+   PUBLIC   dom_nam    ! called by nemogcm::recursive_partition
 
    !! * Control permutation of array indices

branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DOM/domcfg.F90

@@ -86 +86 @@
 !!$      WRITE(*,*) 'ARPDBG: ',narea,': dom_glo: nld{i,j} = ',nldi,nldj
 !!$      WRITE(*,*) 'ARPDBG: ',narea,': dom_glo: nlc{i,j} = ',nlci,nlcj
+!!$      WRITE(*,*) 'ARPDBG: ',narea,': dom_glo: jp{i,j} = ',jpi,jpj
 !!$      WRITE(*,*) 'ARPDBG: ',narea,': dom_glo: n{i,j}mpp = ',nimpp, njmpp
 !!$      WRITE(*,*) 'ARPDBG: ',narea,': dom_glo: jp{i,j}zoom = ',jpizoom, jpjzoom
@@ -110 +111 @@
       !                                   ! (return (m.0,m.1)=(1,0) if data domain gridpoint is to the west/south of the 
       !                                   !local domain, or (m.0,m.1)=(jp.+1,jp.) to the east/north of local domain. 
-#if defined key_mpp_rkpart
+#if 0
+!defined key_mpp_rkpart
       mi0(1:nimpp-1) = 1 !nldi
       DO ji = 0,iesub-1,1

branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

@@ -43 +43 @@
 
    PUBLIC   dom_zgr        ! called by dom_init.F90
+   PUBLIC   zgr_z, zgr_bat, zgr_zco, zgr_zps ! called by nemogcm::recursive_partition
+   PUBLIC   fssig1         ! called by partition_mod::smooth_bathy
 
    !                                       !!* Namelist namzgr_sco *
@@ -54 +56 @@
    !                                        ! ( rn_bb=0; top only, rn_bb =1; top and bottom)
    REAL(wp) ::   rn_hc       =  150._wp     ! Critical depth for s-sigma coordinates
-   PUBLIC rn_sbot_min, rn_sbot_max, rn_theta, rn_thetb, rn_rmax, ln_s_sigma, rn_bb,rn_hc
+   PUBLIC rn_sbot_min, rn_sbot_max, rn_theta, rn_thetb, rn_rmax, &
+          ln_s_sigma, rn_bb, rn_hc
+   PUBLIC ln_zco, ln_zps, ln_sco
+
    !! * Control permutation of array indices
 #  include "oce_ftrans.h90"
@@ -62 +67 @@
 #  include "domzgr_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
+
+   NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco
+   NAMELIST/namzgr_sco/ rn_sbot_max, rn_sbot_min, rn_theta, rn_thetb, &
+                        rn_rmax, ln_s_sigma, rn_bb, rn_hc
+   PUBLIC namzgr, namzgr_sco
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.3.1 , NEMO Consortium (2011)
@@ -88 +98 @@
       INTEGER ::   ioptio = 0   ! temporary integer
       !
-      NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco
-      !!----------------------------------------------------------------------
-
-      REWIND( numnam )                 ! Read Namelist namzgr : vertical coordinate'
+      !NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco
+      !!----------------------------------------------------------------------
+
+      REWIND( numnam )                 ! Read Namelist namzgr : vertical coordinate
       READ  ( numnam, namzgr )
 
@@ -287 +297 @@
 
 
-   SUBROUTINE zgr_bat
+   SUBROUTINE zgr_bat(global_domain)
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE zgr_bat  ***
@@ -317 +327 @@
       !!              - bathy : meter bathymetry (in meters)
       !!----------------------------------------------------------------------
+      LOGICAL, OPTIONAL, INTENT(in) :: global_domain ! Whether dealing with
+                                                     ! whole domain (T) or a 
+                                                     ! sub-domain after domain
+                                                     ! decomposition
+      ! Locals
       INTEGER  ::   ji, jj, jl, jk            ! dummy loop indices
       INTEGER  ::   inum                      ! temporary logical unit
@@ -325 +340 @@
       INTEGER , DIMENSION(jpidta,jpjdta) ::   idta   ! global domain integer data
       REAL(wp), DIMENSION(jpidta,jpjdta) ::   zdta   ! global domain scalar data
+      LOGICAL  ::   is_global 
       !!----------------------------------------------------------------------
 
@@ -330 +346 @@
       IF(lwp) WRITE(numout,*) '    zgr_bat : defines level and meter bathymetry'
       IF(lwp) WRITE(numout,*) '    ~~~~~~~'
+
+      ! Set local flag to signal whether we're dealing with the global domain
+      ! (pre decomposition) or a local part of it. Required by the 
+      ! recursive k-section partitioning.
+      is_global = .FALSE.
+      IF( PRESENT(global_domain) )THEN
+         IF( global_domain ) is_global = .TRUE.
+      END IF
 
       !                                               ! ================== ! 
@@ -347 +371 @@
             ii_bump = jpidta / 2                           ! i-index of the bump center
             ij_bump = jpjdta / 2                           ! j-index of the bump center
-            r_bump  = 50000._wp                            ! bump radius (meters)       
-            h_bump  =  2700._wp                            ! bump height (meters)
+            r_bump  = 0.165*MIN(jpidta,jpjdta)             ! bump radius (grid cells)       
+            h_bump  =  3000._wp                            ! bump height (meters)
             h_oce   = gdepw_0(jpk)                         ! background ocean depth (meters)
             IF(lwp) WRITE(numout,*) '            bump characteristics: '
-            IF(lwp) WRITE(numout,*) '               bump center (i,j)   = ', ii_bump, ii_bump
+            IF(lwp) WRITE(numout,*) '               bump center (i,j)   = ', ii_bump, ij_bump
             IF(lwp) WRITE(numout,*) '               bump height         = ', h_bump , ' meters'
-            IF(lwp) WRITE(numout,*) '               bump radius         = ', r_bump , ' index'
+            IF(lwp) WRITE(numout,*) '               bump radius         = ', r_bump , ' cells'
             IF(lwp) WRITE(numout,*) '            background ocean depth = ', h_oce  , ' meters'
             !                                        
             DO jj = 1, jpjdta                              ! zdta :
                DO ji = 1, jpidta
-                  zi = FLOAT( ji - ii_bump ) * ppe1_m / r_bump
-                  zj = FLOAT( jj - ij_bump ) * ppe2_m / r_bump
+                  !zi = FLOAT( ji - ii_bump ) * ppe1_m / r_bump
+                  !zj = FLOAT( jj - ij_bump ) * ppe2_m / r_bump
+                  zi = FLOAT( ji - ii_bump ) / r_bump
+                  zj = FLOAT( jj - ij_bump ) / r_bump
                   zdta(ji,jj) = h_oce - h_bump * EXP( -( zi*zi + zj*zj ) )
                END DO
@@ -412 +438 @@
          !
          IF( ln_zco )   THEN                          ! zco : read level bathymetry 
-            CALL iom_open ( 'bathy_level.nc', inum )  
-            CALL iom_get  ( inum, jpdom_data, 'Bathy_level', bathy )
+            CALL iom_open ( 'bathy_level.nc', inum ) 
+            IF(is_global)THEN 
+               CALL iom_get ( inum, jpdom_unknown, 'Bathymetry' , bathy, &
+                              kstart=(/jpizoom,jpjzoom/),                &
+                              kcount=(/jpiglo,jpjglo/) )
+            ELSE
+               CALL iom_get  ( inum, jpdom_data, 'Bathy_level', bathy )
+            END IF
+
             CALL iom_close( inum )
             mbathy(:,:) = INT( bathy(:,:) )
@@ -446 +479 @@
          IF( ln_zps .OR. ln_sco )   THEN              ! zps or sco : read meter bathymetry
             CALL iom_open ( 'bathy_meter.nc', inum ) 
-            CALL iom_get  ( inum, jpdom_data, 'Bathymetry', bathy )
+            IF(is_global)THEN 
+               CALL iom_get ( inum, jpdom_unknown, 'Bathymetry' , bathy, &
+                              kstart=(/jpizoom,jpjzoom/),                &
+                              kcount=(/jpiglo,jpjglo/) )
+            ELSE
+               CALL iom_get  ( inum, jpdom_data, 'Bathymetry', bathy )
+            END IF
             CALL iom_close( inum )
             !                                                ! =====================
@@ -516 +555 @@
          zhmin = gdepw_0(ik+1)                                                         ! minimum depth = ik+1 w-levels 
          WHERE( bathy(:,:) <= 0._wp )   ;   bathy(:,:) = 0._wp                         ! min=0     over the lands
-         ELSE WHERE                     ;   bathy(:,:) = MAX(  zhmin , bathy(:,:)  )   ! min=zhmin over the oceans
+         ELSEWHERE                      ;   bathy(:,:) = MAX(  zhmin , bathy(:,:)  )   ! min=zhmin over the oceans
          END WHERE
          IF(lwp) write(numout,*) 'Minimum ocean depth: ', zhmin, ' minimum number of ocean levels : ', ik
@@ -730 +769 @@
       !!                                     (min value = 1 over land)
       !!----------------------------------------------------------------------
+      !USE arpdebugging, ONLY: dump_array
       USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
       USE wrk_nemo, ONLY:   zmbk => wrk_2d_1
@@ -744 +784 @@
       IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~'
       !
+      !CALL dump_array(0, 'mbathy', mbathy, withHalos=.TRUE.)
+
       mbkt(:,:) = MAX( mbathy(:,:) , 1 )    ! bottom k-index of T-level (=1 over land)
       !                                     ! bottom k-index of W-level = mbkt+1
@@ -755 +797 @@
       zmbk(:,:) = REAL( mbku(:,:), wp )   ;   CALL lbc_lnk(zmbk,'U',1.)   ;   mbku  (:,:) = MAX( INT( zmbk(:,:) ), 1 )
       zmbk(:,:) = REAL( mbkv(:,:), wp )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv  (:,:) = MAX( INT( zmbk(:,:) ), 1 )
+      !
+      ! Compute and store the deepest bottom level of any grid-type at each grid point
+      ! For use in removing data below ocean floor from halo exchanges.
+      mbkmax(:,:) = MAX(mbkt(:,:)+1, mbku(:,:), mbkv(:,:))
       !
       IF( wrk_not_released(2, 1) )   CALL ctl_stop('zgr_bot_level: failed to release workspace array')
@@ -809 +855 @@
 
 
-   SUBROUTINE zgr_zps
+   SUBROUTINE zgr_zps(pre_domain_decomp)
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE zgr_zps  ***
@@ -858 +904 @@
       !! DCSE_NEMO: wrk_3d_1 renamed, need additional directive
 !FTRANS zprt :I :I :z
+      LOGICAL, INTENT(in), OPTIONAL :: pre_domain_decomp
       !!
       INTEGER  ::   ji, jj, jk       ! dummy loop indices
@@ -892 +939 @@
       bathy(:,:) = MIN( zmax ,  bathy(:,:) )    ! bounded value of bathy (min already set at the end of zgr_bat)
       WHERE( bathy(:,:) == 0._wp )   ;   mbathy(:,:) = 0       ! land  : set mbathy to 0
-      ELSE WHERE                     ;   mbathy(:,:) = jpkm1   ! ocean : initialize mbathy to the max ocean level
+      ELSEWHERE                     ;   mbathy(:,:) = jpkm1   ! ocean : initialize mbathy to the max ocean level
       END WHERE
 
@@ -903 +950 @@
          WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth )   mbathy(:,:) = jk-1
       END DO
+
+      ! If we've been called before domain decomposition then we only want to compute
+      ! mbathy and the return.
+      IF( PRESENT(pre_domain_decomp) )THEN
+         IF( pre_domain_decomp )RETURN
+      ENDIF
 
       ! Scale factors and depth at T- and W-points
@@ -1208 +1261 @@
       USE mapcomm_mod, ONLY: trimmed, cyclic_bc
       USE mapcomm_mod, ONLY: nidx, eidx, sidx, widx
-!      USE arpdebugging, ONLY: dump_array
+
       !! DCSE_NEMO: wrk_nemo module variables renamed, need additional directives
 !FTRANS gsigw3 :I :I :z
@@ -1227 +1280 @@
       !
 
-      NAMELIST/namzgr_sco/ rn_sbot_max, rn_sbot_min, rn_theta, rn_thetb, rn_rmax, ln_s_sigma, rn_bb, rn_hc
+!      NAMELIST/namzgr_sco/ rn_sbot_max, rn_sbot_min, rn_theta, rn_thetb, rn_rmax, ln_s_sigma, rn_bb, rn_hc
       !!----------------------------------------------------------------------
 
@@ -1319 +1372 @@
                                                          jl, zrmax, INT( SUM(zmsk(:,:) ) )
          !
-!!$         IF(jl < 6)THEN ! .OR. (MOD(jl,1000) == 0) )THEN
-!!$            CALL dump_array(jl, 'zenv_before', zenv, withHalos=.TRUE.)
-!!$            CALL dump_array(jl, 'ztmp_before', ztmp, withHalos=.TRUE.)
-!!$            CALL dump_array(jl, 'zmsk_before', zmsk, withHalos=.TRUE.)
-!!$         END IF
 
          ! Copy current surface before next smoothing iteration 
@@ -1361 +1409 @@
          ! Apply lateral boundary condition but do not zero on closed boundaries
          CALL lbc_lnk( zenv, 'T', 1._wp, lzero=.FALSE. )
-
-!!$         IF(jl < 6)THEN ! .OR. (MOD(jl,1000) == 0) )THEN
-!!$            CALL dump_array(jl, 'zenv', zenv, withHalos=.TRUE.)
-!!$            CALL dump_array(jl, 'ztmp', ztmp, withHalos=.TRUE.)
-!!$            CALL dump_array(jl, 'zmsk', zmsk, withHalos=.TRUE.)
-!!$         END IF
 
          !                                                  ! ================ !
@@ -1709 +1751 @@
                   CALL ctl_stop( ctmp1 )
                ENDIF
+#if defined key_vvl
                IF( gdepw_1(ji,jj,jk) < 0._wp .OR. gdept_1(ji,jj,jk) < 0._wp ) THEN
                   WRITE(ctmp1,*) 'zgr_sco :   gdepw or gdept =< 0  at point (i,j,k)= ', ji, jj, jk
                   CALL ctl_stop( ctmp1 )
                ENDIF
+#endif
             END DO
          END DO

branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90

@@ -166 +166 @@
          ENDIF                                                        ! explicit case not coded yet with AGRIF
       ENDIF
-      !
+
    END SUBROUTINE istate_init
 

branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90

@@ -118 +118 @@
          WRITE(numout,*)      '           ','   jpij    : ', jpij
          WRITE(numout,*) '          mpp local domain info (mpp)'
+#if defined key_mpp_rkpart
+         WRITE(numout,*) '             recursive k-section decomposition used. See file: domain_decomp.ps'
+         WRITE(numout,*) '             jpreci  : ', jpreci, '   jprecj  : ', jprecj
+#else
+         ! These parameters are only used in the original, regular domain
+         ! decomposition scheme so we don't print them if we're using
+         ! recursive k-section partitioning.
          WRITE(numout,*) '             jpni    : ', jpni, '   jpreci  : ', jpreci
          WRITE(numout,*) '             jpnj    : ', jpnj, '   jprecj  : ', jprecj
          WRITE(numout,*) '             jpnij   : ', jpnij
+#endif
          WRITE(numout,*) '          lateral domain boundary condition type : jperio  = ', jperio
          WRITE(numout,*)

branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90

@@ -102 +102 @@
       USE oce     , ONLY:   ze3u_f => ta       , ze3v_f => sa       ! (ta,sa) used as 3D workspace
       USE wrk_nemo, ONLY:   zs_t   => wrk_2d_1 , zs_u_1 => wrk_2d_2 , zs_v_1 => wrk_2d_3
+      USE arpdebugging, ONLY: dump_array
       !! DCSE_NEMO: need additional directives for renamed module variables
 !FTRANS ze3u_f :I :I :z
@@ -126 +127 @@
          IF(lwp) WRITE(numout,*) '~~~~~~~'
       ENDIF
+
+!      CALL dump_array(kt, 'ua_nxt_start',ua(:,:,1),withHalos=.TRUE.)
 
 #if defined key_dynspg_flt

branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90

@@ -133 +133 @@
          END DO
       ENDIF
-
 
       SELECT CASE ( nspg )                       ! compute surf. pressure gradient trend and add it to the general trend

branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_flt.F90

@@ -112 +112 @@
       !!---------------------------------------------------------------------
       USE oce, ONLY:   zub   => ta , zvb   => sa   ! (ta,sa) used as workspace
+!      USE arpdebugging, ONLY: dump_array
       !! DCSE_NEMO: need additional directives for renamed module variables
 !FTRANS zub :I :I :z
@@ -137 +138 @@
          !                                                        ! gcx, gcxb
       ENDIF
+
+!      CALL dump_array(kt, 'spgu',spgu,withHalos=.TRUE.)
+!      CALL dump_array(kt, 'sshn',sshn,withHalos=.TRUE.)
+!#if defined key_z_first
+!      CALL dump_array(kt, 'ua',ua(1,:,:),withHalos=.TRUE.)
+!#else
+!      CALL dump_array(kt, 'ua',ua(:,:,1),withHalos=.TRUE.)
+!#endif
 
       ! Local constant initialization

branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90

@@ -157 +157 @@
          CALL trd_mod(ztrdu, ztrdv, jpdyn_trd_zad, 'DYN', kt)
       ENDIF
+
       !                             ! Control print
       IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' zad  - Ua: ', mask1=umask,   &

branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90

@@ -63 +63 @@
       USE oce     , ONLY:  zwd  => ta       , zws   => sa   ! (ta,sa) used as 3D workspace
       USE wrk_nemo, ONLY:   zwi => wrk_3d_3                 ! 3D workspace
+      USE arpdebugging, ONLY: dump_array
       !! DCSE_NEMO: need additional directives for renamed module variables
 !FTRANS zwd :I :I :z
@@ -88 +89 @@
       ! --------------------------------
       z1_p2dt = 1._wp / p2dt      ! inverse of the timestep
+
+      !CALL dump_array(kt, 'utau_pre_zdf',utau(:,:),withHalos=.TRUE.)
+      !CALL dump_array(kt, 'utaub_pre_zdf',utau_b(:,:),withHalos=.TRUE.)
 
       ! 1. Vertical diffusion on u

branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/LBC/exchmod.F90

@@ -11 +11 @@
   ! Make some key parameters from mapcomm_mod available to all who 
   ! USE this module
-  USE mapcomm_mod, ONLY: Iminus, Iplus, Jminus, Jplus, NONE
+  USE mapcomm_mod, ONLY: Iminus, Iplus, Jminus, Jplus, NONE, &
+                         jeub
   IMPLICIT none
 
@@ -63 +64 @@
   INTEGER, DIMENSION(:), ALLOCATABLE, SAVE :: nrank_north  ! dim. ndim_rank_north, number
                                                            ! of the procs belonging to ncomm_north
+  LOGICAL, SAVE :: do_nfold ! Whether this PE contributes to N-fold exchange
+                            !  -  takes domain trimming into account.
   INTEGER, PARAMETER :: num_nfold_rows = 4 ! No. of rows at the top of the 
                                            ! global domain to use in applying 
                                            ! the north-fold condition (no value
-                                           ! other than 4 currently supported)
+                                           ! other than 4 currently tested)
+
+  INTEGER, DIMENSION(:), ALLOCATABLE, SAVE :: nfold_npts ! How many points each
+                                                         ! northern proc contrib
+                                                         ! to nfold exchange
 
 !FTRANS r3dptr :I :I :z
@@ -112 +119 @@
          nrank_north, north_root, ndim_rank_north, &
          ngrp_north, ngrp_world, ncomm_north, &
+         num_nfold_rows, do_nfold, nfold_npts, &
          exchmod_alloc, add_exch, bound_exch_list, &
-         Iminus, Iplus, Jminus, Jplus, NONE, num_nfold_rows, &
-         lbc_exch3, lbc_exch2, & !lbc_exch3i, lbc_exch2i, &
-         MPI_COMM_WORLD, MPI_Wtime
+         Iminus, Iplus, Jminus, Jplus, NONE, &
+         lbc_exch3, lbc_exch2
+
+#if defined key_mpp_mpi
+  PUBLIC MPI_COMM_WORLD, MPI_Wtime
+#endif
 
   ! MPI only
@@ -285 +296 @@
 
   SUBROUTINE bound_exch_generic ( b2, ib2, b3, ib3, nhalo, nhexch, &
-       comm1, comm2, comm3, comm4,      &
-       cd_type, lfill, isgn, lzero )
+                                  comm1, comm2, comm3, comm4,      &
+                                  cd_type, lfill, pval, isgn, lzero )
     USE par_oce, ONLY: wp, jpreci, jprecj, jpim1
     USE dom_oce, ONLY: nlci, nlcj, nldi, nlei, nldj, nlej, &
-         nperio, nbondi, npolj
+                       nperio, nbondi, npolj, narea
     USE mapcomm_mod, ONLY: Iminus, Iplus, NONE, ilbext, iubext, cyclic_bc
+    USE mapcomm_mod, ONLY: trimmed, eidx, widx
     IMPLICIT none
     INTEGER, INTENT(in)  :: nhalo,nhexch
@@ -302 +314 @@
     CHARACTER(len=1),  INTENT(in) :: cd_type
     LOGICAL, OPTIONAL, INTENT(in) :: lfill
+    REAL(wp),OPTIONAL, INTENT(in) :: pval  ! background value (used at closed boundaries)
     INTEGER, OPTIONAL, INTENT(in) :: isgn
-    LOGICAL, OPTIONAL, INTENT(in) :: lzero ! Whether to zero halos on closed boundaries
+    LOGICAL, OPTIONAL, INTENT(in) :: lzero ! Whether to set halo values on closed boundaries
     ! Local arguments
     INTEGER :: itag          ! Communication handle
@@ -310 +323 @@
     INTEGER :: ileft, iright ! First and last x-coord of internal points
     INTEGER :: kdim1
+    INTEGER  :: iland ! Land values - zero by default unless pval passed in.
+    REAL(wp) :: zland !  "     "
     LOGICAL :: lfillarg, lzeroarg
     !!--------------------------------------------------------------------
@@ -324 +339 @@
     lfillarg = .FALSE.
     isgnarg = 1
+    zland = 0.0_wp
 
     IF( PRESENT(lfill) ) lfillarg = lfill
     IF( PRESENT(isgn)  ) isgnarg  = isgn
     IF( PRESENT(lzero) ) lzeroarg = lzero
+    IF( PRESENT(pval)  ) zland    = pval 
+    iland=INT(zland)
 
     ! Find out the size of 3rd dimension of the array
@@ -356 +374 @@
        ! have cyclic E-W boundary conditions.
        ileft = nldi
-       IF(ilbext .AND. cyclic_bc)ileft = ileft + 1
+       IF( (ilbext .AND. (.NOT. trimmed(widx,narea))) .AND. cyclic_bc) &
+                                                     ileft = ileft + 1
 
        iright = nlei
-       IF(iubext .AND. cyclic_bc)iright = iright - 1
+       IF( (iubext .AND. (.NOT. trimmed(eidx,narea))) .AND. cyclic_bc) &
+                                                    iright = iright - 1
 
        IF ( PRESENT(b2) ) THEN
@@ -527 +547 @@
        END IF
 
-    ELSE ! lfillarg is .FALSE.
+    ELSE ! lfillarg is .FALSE. - standard closed or cyclic treatment
 
        !                                        ! East-West boundaries
        !                                        ! ====================
+       !   nbondi == 2 when a single sub-domain spans the whole width
+       !   of the global domain
        IF( nbondi == 2 .AND.   &      ! Cyclic east-west
             &    (nperio == 1 .OR. nperio == 4 .OR. nperio == 6) ) THEN
@@ -568 +590 @@
                 SELECT CASE ( cd_type )
                 CASE ( 'T', 'U', 'V', 'W' , 'I' )
-                   b2(1:jpreci         , :) = 0._wp ! Western halo
-                   b2(nlci-jpreci+1:jpi, :) = 0._wp ! Eastern halo
+                   b2(1:jpreci         , :) = zland ! Western halo
+                   b2(nlci-jpreci+1:jpi, :) = zland ! Eastern halo
                 CASE ( 'F' )
-                   b2(nlci-jpreci+1:jpi, :) = 0._wp ! Eastern halo
+                   b2(nlci-jpreci+1:jpi, :) = zland ! Eastern halo
                 END SELECT
              ELSE IF ( PRESENT(ib2) ) THEN
                 SELECT CASE ( cd_type )
                 CASE ( 'T', 'U', 'V', 'W' , 'I' )
-                   ib2(1:jpreci         , :) = 0 ! Western halo
-                   ib2(nlci-jpreci+1:jpi, :) = 0 ! Eastern halo
+                   ib2(1:jpreci         , :) = iland ! Western halo
+                   ib2(nlci-jpreci+1:jpi, :) = iland ! Eastern halo
                 CASE ( 'F' )
-                   ib2(nlci-jpreci+1:jpi, :) = 0 ! Eastern halo
+                   ib2(nlci-jpreci+1:jpi, :) = iland ! Eastern halo
                 END SELECT
              ELSE IF ( PRESENT(b3) ) THEN
@@ -588 +610 @@
                       DO ji=1,jpreci,1
                          DO jk=1,jpk,1
-                            b3(ji, jj, jk) = 0._wp
+                            b3(ji, jj, jk) = zland
                          END DO
                       END DO
                       DO ji=nlci-jpreci+1,jpi,1
                          DO jk=1,jpk,1
-                            b3(ji, jj, jk) = 0._wp
+                            b3(ji, jj, jk) = zland
                          END DO
                       END DO
                    END DO
 #else
-                   b3(1:jpreci         , :, :) = 0._wp
-                   b3(nlci-jpreci+1:jpi, :, :) = 0._wp
+                   b3(1:jpreci         , :, :) = zland
+                   b3(nlci-jpreci+1:jpi, :, :) = zland
 #endif
                 CASE ( 'F' )
@@ -606 +628 @@
                       DO ji = nlci-jpreci+1,jpi,1
                          DO jk = 1,jpk,1
-                            b3(ji, jj, jk) = 0._wp
+                            b3(ji, jj, jk) = zland
                          END DO
                       END DO
                    END DO
 #else
-                   b3(nlci-jpreci+1:jpi, :, :) = 0._wp
+                   b3(nlci-jpreci+1:jpi, :, :) = zland
 #endif
                 END SELECT
@@ -617 +639 @@
                 SELECT CASE ( cd_type )
                 CASE ( 'T', 'U', 'V', 'W' )
-                   ib3(1:jpreci         , :, :) = 0
-                   ib3(nlci-jpreci+1:jpi, :, :) = 0
+                   ib3(1:jpreci         , :, :) = iland
+                   ib3(nlci-jpreci+1:jpi, :, :) = iland
                 CASE ( 'F' )
-                   ib3(nlci-jpreci+1:jpi, :, :) = 0
+                   ib3(nlci-jpreci+1:jpi, :, :) = iland
                 END SELECT
              END IF
@@ -630 +652 @@
        IF( lzeroarg )THEN
 
-          !                                        ! North-South boundaries
-          !                                        ! ======================
+          !                             ! North-South boundaries (always closed)
+          !                             ! ======================
           IF ( PRESENT(b2) ) THEN
              SELECT CASE ( cd_type )
              CASE ( 'T', 'U', 'V', 'W' , 'I' )
-                b2(:,1:nldj-1         ) = 0._wp
-                b2(:,nlcj-jprecj+1:jpj) = 0._wp
+                !b2(:,1:nldj-1         ) = zland
+                ! Below is what is done in original lib_mpp.F90
+                b2(:,1:jprecj         ) = zland
+                b2(:,nlcj-jprecj+1:jpj) = zland
              CASE ( 'F' )
-                b2(:,nlcj-jprecj+1:jpj) = 0._wp
+                b2(:,nlcj-jprecj+1:jpj) = zland
              END SELECT
           ELSE IF ( PRESENT(ib2) ) THEN
              SELECT CASE ( cd_type )
              CASE ( 'T', 'U', 'V', 'W' , 'I' )
-                ib2(:,1:jprecj         ) = 0
-                ib2(:,nlcj-jprecj+1:jpj) = 0
+                ib2(:,1:jprecj         ) = iland
+                ib2(:,nlcj-jprecj+1:jpj) = iland
              CASE ( 'F' )
-                ib2(:,nlcj-jprecj+1:jpj) = 0
+                ib2(:,nlcj-jprecj+1:jpj) = iland
              END SELECT
           ELSE IF ( PRESENT(b3) ) THEN
@@ -652 +676 @@
              CASE ( 'T', 'U', 'V', 'W' )
 #if defined key_z_first
-                DO jj=1,nldj-1,1
+                DO jj=1,jprecj,1
                    DO ji=1,jpi,1
                       DO jk = 1,jpk,1
-                         b3(ji, jj, jk) = 0._wp
+                         b3(ji, jj, jk) = zland
                       END DO
                    END DO
@@ -662 +686 @@
                    DO ji=1,jpi,1
                       DO jk = 1,jpk,1
-                         b3(ji, jj, jk) = 0._wp
+                         b3(ji, jj, jk) = zland
                       END DO
                    END DO
                 END DO
 #else
-                b3(:, 1:nldj-1         , :) = 0._wp
-                b3(:, nlcj-jprecj+1:jpj, :) = 0._wp
+                b3(:, 1:jprecj         , :) = zland
+                b3(:, nlcj-jprecj+1:jpj, :) = zland
 #endif
              CASE ( 'F' )
@@ -675 +699 @@
                    DO ji=1,jpi,1
                       DO jk = 1,jpk,1
-                         b3(ji, jj, jk) = 0._wp
+                         b3(ji, jj, jk) = zland
                       END DO
                    END DO
                 END DO
 #else
-                b3(:, nlcj-jprecj+1:jpj, :) = 0._wp
+                b3(:, nlcj-jprecj+1:jpj, :) = zland
 #endif
              END SELECT
@@ -686 +710 @@
              SELECT CASE ( cd_type )
              CASE ( 'T', 'U', 'V', 'W' )
-                ib3(:, 1:jprecj         , :) = 0
-                ib3(:, nlcj-jprecj+1:jpj, :) = 0
+                ib3(:, 1:jprecj         , :) = iland
+                ib3(:, nlcj-jprecj+1:jpj, :) = iland
              CASE ( 'F' )
-                ib3(:, nlcj-jprecj+1:jpj, :) = 0
+                ib3(:, nlcj-jprecj+1:jpj, :) = iland
              END SELECT
           END IF
@@ -726 +750 @@
        ! We only need to repeat the East and West halo swap if there
        ! IS a north-fold in the configuration.
-       SELECT CASE (npolj)
-
-       CASE ( 3, 4, 5, 6 )
-
-          ! Update East and West halos as required
+       !SELECT CASE (npolj)
+
+       !CASE ( 3, 4, 5, 6 )
+       IF(ndim_rank_north > 0)THEN
+
+          ! Update East and West halos as required - no data sent north
+          ! as it's only the northern-most PEs that have been affected 
+          ! by the north-fold condition.
           ! ARPDBG - inefficient since all PEs do halo swap and only 
           ! those affected by the north fold actually need to - can 
           ! this be done within apply_north_fold?
           CALL exchs_generic (b2=b2,ib2=ib2,b3=b3,ib3=ib3, nhalo=nhalo, &
-               nhexch=nhexch, handle=itag,               &
-               comm1=Iplus,comm2=Iminus,comm3=NONE,comm4=NONE, &
-               cd_type=cd_type, lfill=lfillarg)
+                              nhexch=nhexch, handle=itag,               &
+                              comm1=Iplus,comm2=Iminus,comm3=Jplus,comm4=Jminus, &
+                              cd_type=cd_type, lfill=lfillarg)
 
           !CALL exchr_generic (b2=b2,ib2=ib2,b3=b3,ib3=ib3,nhalo=nhalo, &
@@ -743 +770 @@
           !                    comm1=Iplus,comm2=Iminus,comm3=NONE,comm4=NONE)
           !                           comm1=comm1,comm2=comm2,comm3=comm3,comm4=comm4 )
-       END SELECT    ! npolj 
+       END IF        ! ndim_rank_north > 0
+       !END SELECT    ! npolj 
 
     END IF
@@ -1160 +1188 @@
 
        DO ifield = 1, nfields, 1
-          IF( npolj /= 0 )THEN ! only for northern procs.
+          IF( npolj /= 0 .AND. do_nfold )THEN ! only for northern procs.
 
              IF(ASSOCIATED(list(ifield)%r2dptr))THEN
@@ -1180 +1208 @@
        END DO
 
-!!$       IF( npolj /= 0 ) CALL mpp_lbc_north_list( list, nfields ) ! only for northern procs.
+!!$       IF( npolj /= 0 .AND. do_nfold ) CALL mpp_lbc_north_list( list, nfields ) ! only for northern procs.
 
     END SELECT   ! jpni
@@ -1943 +1971 @@
 
     CASE DEFAULT   ! more than 1 proc along I
-       IF( npolj /= 0 )   CALL mpp_lbc_north( b2, cd_type, psgn )   ! only for northern procs.
+       IF( npolj /= 0 .AND. do_nfold ) CALL mpp_lbc_north( b2, cd_type, psgn )   ! only for northern procs.
 
     END SELECT   ! jpni
@@ -2074 +2102 @@
 
     CASE DEFAULT   ! more than 1 proc along I
-       IF( npolj /= 0 )   CALL mpp_lbc_north( ib2, cd_type, isgn )   ! only for northern procs.
+       IF( npolj /= 0 .AND. do_nfold) CALL mpp_lbc_north( ib2, cd_type, isgn )   ! only for northern procs.
 
     END SELECT   ! jpni
@@ -2285 +2313 @@
 
     CASE DEFAULT ! more than 1 proc along I
-       IF ( npolj /= 0 ) CALL mpp_lbc_north (b3, cd_type, psgn) ! only for northern procs.
+       IF ( npolj /= 0 .AND. do_nfold) CALL mpp_lbc_north (b3, cd_type, psgn) ! only for northern procs.
 
     END SELECT ! jpni 
@@ -2493 +2521 @@
 
     CASE DEFAULT ! more than 1 proc along I
-       IF ( npolj /= 0 ) CALL mpp_lbc_north ( ib3, cd_type, isgn) ! only for northern procs.
+       IF ( npolj /= 0 .AND. do_nfold) CALL mpp_lbc_north ( ib3, cd_type, isgn) ! only for northern procs.
 
     END SELECT ! jpni 
@@ -2557 +2585 @@
     ELSE
        ! This section is both for error checking and allows me to be lazy in the 
-       ! testing code - I don't have to check which arrays I've been passed.
+       ! testing code - I don't have to check which arrays I've been passed 
+       ! before I call this routine.
        WRITE (*,*) 'WARNING: add_exch called without a ptr to an array - will be ignored'
        RETURN
@@ -2588 +2617 @@
    SUBROUTINE bound_exch2 (b, nhalo, nhexch,           &
                            comm1, comm2, comm3, comm4, &
-                           cd_type, lfill, isgn, lzero )
+                           cd_type, lfill, pval, isgn, lzero )
       !!----------------------------------------------------------------------
       !!----------------------------------------------------------------------
@@ -2600 +2629 @@
       INTEGER, OPTIONAL, INTENT(in) :: isgn
       LOGICAL, OPTIONAL, INTENT(in) :: lzero
+      REAL(wp),OPTIONAL, INTENT(in) :: pval
 
       CALL bound_exch_generic( b2=b,nhalo=nhalo,nhexch=nhexch, &
               comm1=comm1,comm2=comm2,comm3=comm3,comm4=comm4, &
-              cd_type=cd_type, lfill=lfill, isgn=isgn, lzero=lzero )
+              cd_type=cd_type, lfill=lfill, pval=pval,         &
+              isgn=isgn, lzero=lzero )
       RETURN
    END SUBROUTINE bound_exch2
@@ -2609 +2640 @@
 
    SUBROUTINE bound_exch2i (b, nhalo, nhexch, comm1, comm2, comm3, comm4, &
-                            cd_type, lfill, isgn, lzero )
+                            cd_type, lfill, pval, isgn, lzero )
       !!----------------------------------------------------------------------
       !!----------------------------------------------------------------------
@@ -2621 +2652 @@
       INTEGER, OPTIONAL, INTENT(in) :: isgn
       LOGICAL, OPTIONAL, INTENT(in) :: lzero
+      REAL(wp),OPTIONAL, INTENT(in) :: pval
 
       CALL bound_exch_generic (ib2=b,nhalo=nhalo,nhexch=nhexch,           &
                          comm1=comm1,comm2=comm2,comm3=comm3,comm4=comm4, &
-                         cd_type=cd_type, lfill=lfill, isgn=isgn, lzero=lzero )
+                         cd_type=cd_type, lfill=lfill, pval=pval,         &
+                         isgn=isgn, lzero=lzero )
       RETURN
    END SUBROUTINE bound_exch2i
@@ -2630 +2663 @@
 
    SUBROUTINE bound_exch3 (b, nhalo, nhexch, comm1, comm2, comm3, &
-                          comm4, cd_type, lfill, isgn, lzero)
+                          comm4, cd_type, lfill, pval, isgn, lzero)
       !!----------------------------------------------------------------------
       !!----------------------------------------------------------------------
@@ -2642 +2675 @@
       INTEGER, OPTIONAL, INTENT(in) :: isgn
       LOGICAL, OPTIONAL, INTENT(in) :: lzero
+      REAL(wp),OPTIONAL, INTENT(in) :: pval
 
       CALL bound_exch_generic ( b3=b,nhalo=nhalo,nhexch=nhexch,&
               comm1=comm1,comm2=comm2,comm3=comm3,comm4=comm4, &
-              cd_type=cd_type, lfill=lfill, isgn=isgn, lzero=lzero )
+              cd_type=cd_type, lfill=lfill, pval=pval,         &
+              isgn=isgn, lzero=lzero )
       RETURN
    END SUBROUTINE bound_exch3
@@ -2651 +2686 @@
 
    SUBROUTINE bound_exch3i (b, nhalo, nhexch, comm1, comm2, comm3, &
-                           comm4, cd_type, lfill, isgn, lzero)
+                           comm4, cd_type, lfill, pval, isgn, lzero)
       !!----------------------------------------------------------------------
       !!----------------------------------------------------------------------
@@ -2662 +2697 @@
       INTEGER, OPTIONAL, INTENT(in) :: isgn
       LOGICAL, OPTIONAL, INTENT(in) :: lzero
+      REAL(wp),OPTIONAL, INTENT(in) :: pval
 
       CALL bound_exch_generic ( ib3=b,nhalo=nhalo,nhexch=nhexch, &
                 comm1=comm1,comm2=comm2,comm3=comm3,comm4=comm4, &
-                cd_type=cd_type, lfill=lfill, isgn=isgn, lzero=lzero )
+                cd_type=cd_type, lfill=lfill, pval=pval,         &
+                isgn=isgn, lzero=lzero )
 
    END SUBROUTINE bound_exch3i
@@ -2695 +2732 @@
       LOGICAL :: lfill
 
-      ! ARPDBG - don't know whether pval currently maps into exchmod framework
-      IF(PRESENT(pval))THEN
-         CALL ctl_stop('STOP','lbc_exch2: got pval argument - NOT IMPLEMENTED')
-         RETURN
-      END IF
-
       lfill = .FALSE.
       IF(PRESENT(cd_mpp))THEN
@@ -2708 +2739 @@
       CALL bound_exch_generic( b2=pt2d,nhalo=jpreci,nhexch=jpreci, &
             comm1=Iplus,comm2=Iminus,comm3=Jplus,comm4=Jminus, &
-            cd_type=cd_type, lfill=lfill, isgn=INT(psgn), lzero=lzero )
+            cd_type=cd_type, lfill=lfill, pval=pval, isgn=INT(psgn), &
+            lzero=lzero )
 
    END SUBROUTINE lbc_exch2
@@ -2731 +2763 @@
       LOGICAL :: lfill
 
-      ! ARPDBG - don't know whether pval currently maps into exchmod framework
-      IF(PRESENT(pval))THEN
-         CALL ctl_stop('STOP','lbc_exch3: got pval argument - NOT IMPLEMENTED')
-         RETURN
-      END IF
-
       lfill = .FALSE.
       IF(PRESENT(cd_mpp))THEN
@@ -2742 +2768 @@
       END IF
 
-      CALL bound_exch_generic ( b3=ptab3d,nhalo=jpreci,nhexch=jpreci,&
-             comm1=Iplus,comm2=Iminus,comm3=Jplus,comm4=Jminus, &
-             cd_type=cd_type, lfill=lfill, isgn=INT(psgn), lzero=lzero )
+      CALL bound_exch_generic ( b3=ptab3d,nhalo=jpreci,nhexch=jpreci, &
+             comm1=Iplus,comm2=Iminus,comm3=Jplus,comm4=Jminus,       &
+             cd_type=cd_type, lfill=lfill, pval=pval, isgn=INT(psgn), &
+             lzero=lzero )
 
    END SUBROUTINE lbc_exch3
@@ -2773 +2800 @@
     USE mapcomm_mod, ONLY: Iplus, Jplus, Iminus, Jminus, IplusJplus,       &
                            IminusJminus, IplusJminus, IminusJplus,  &
-                           nsend, nxsend, nysend, nxsendp,nysendp,nsendp, &
+                           nsend, nxsend, nysend, nxsendp,nysendp,nzsendp, &
+                           nsendp, &
                            destination,dirsend, dirrecv,                  &
                            isrcsendp,jsrcsendp, idesrecvp, jdesrecvp,     &
-                           nrecv, nxrecv,nyrecv,nxrecvp,nyrecvp,nrecvp,   &
+                           nrecv,  &
+                           nxrecvp,nyrecvp,nzrecvp, nrecvp, nrecvp2d,  &
                            source, iesub, jesub,  &
                            MaxCommDir, MaxComm, cyclic_bc,      &
                            nrecvp, npatchsend, npatchrecv
-    USE lib_mpp,     ONLY: mpi_comm_opa, ctl_stop
+    USE lib_mpp,     ONLY: ctl_stop
+#if defined key_mpp_mpi
+    USE lib_mpp,     ONLY: mpi_comm_opa
+#endif
 #if ( defined DEBUG && defined DEBUG_EXCHANGE ) || defined DEBUG_COMMS
     USE dom_oce,     ONLY: narea
@@ -2793 +2825 @@
 
     LOGICAL :: enabled(0:MaxCommDir, maxExchItems)
-    INTEGER :: ides, ierr, irecv, isend, &
-               isrc, jdes, jsrc, nxr, nyr,        &
-               nxs, nys, tag, tag_orig,           &
+    INTEGER :: ides, ierr, irecv, isend,        &
+               isrc, jdes, jsrc, tag, tag_orig, &
                ibeg, iend, jbeg, jend
     INTEGER :: i, j, k, ic, ifield, ipatch ! Loop counters
@@ -2801 +2832 @@
     INTEGER :: npacked
     INTEGER :: handle
+#if defined key_mpp_mpi
     INTEGER :: status(MPI_status_size)
     INTEGER :: astatus(MPI_status_size,MaxComm)
+#endif
     INTEGER :: r2dcount, r3dcount, i2dcount, i3dcount
     ! Indices into int and real copy buffers
@@ -2826 +2859 @@
 #endif
 
-    CALL prof_region_begin(ARPEXCHS_LIST, "Exchs_list", iprofStat)
+    !CALL prof_region_begin(ARPEXCHS_LIST, "Exchs_list", iprofStat)
 
     ! Allocate a communications tag/handle and a flags array.
@@ -2841 +2874 @@
        ! Check halo width is in range.
        IF ( list(ifield)%halo_width.GT.jpreci ) THEN
-          CALL ctl_stop('STOP','exchs: halo width greater than maximum')
+          CALL ctl_stop('STOP', &
+                        'exchs_generic_list: halo width greater than maximum')
           RETURN
        ENDIF
@@ -2881 +2915 @@
     IF( have_real_field )THEN
 
-       ALLOCATE(recvBuff(jpkdta*maxrecvpts*nfields,nrecv),stat=ierr)
+       ALLOCATE(recvBuff(maxrecvpts*nfields,nrecv),stat=ierr)
        !WRITE(*,"('Allocated ',I7,' reals for recv buff')") &
        !                                 jpkdta*maxrecvpts*nfields
@@ -2898 +2932 @@
     IF( have_int_field .AND. (ierr == 0) )THEN
 
-       ALLOCATE(recvIBuff(jpkdta*maxrecvpts*nfields,nrecv),stat=ierr)
+       ALLOCATE(recvIBuff(maxrecvpts*nfields,nrecv),stat=ierr)
        !WRITE(*,"('Allocated ',I7,' ints for recv buff')") &
        !                                 jpkdta*maxrecvpts*nfields
@@ -2927 +2961 @@
        i3dcount = 0
 
-       IF(source(irecv).GE.0 .AND. nrecvp(irecv,1).GT.0 ) THEN
+       IF( source(irecv).GE.0 .AND. &
+           ( (nrecvp(irecv,1) > 0) .OR. (nrecvp2d(irecv,1) > 0) ) ) THEN
 
           ! This loop is to allow for different fields to have different
@@ -2935 +2970 @@
              IF ( enabled(dirrecv(irecv), ifield) ) THEN
                 IF( ASSOCIATED(list(ifield)%r2dptr) )THEN
-                   r2dcount = r2dcount + 1
+                   r2dcount = r2dcount + nrecvp2d(irecv,1)
                 ELSE IF( ASSOCIATED(list(ifield)%i2dptr) )THEN
-                   i2dcount = i2dcount + 1
+                   i2dcount = i2dcount + nrecvp2d(irecv,1)
                 ELSE IF( ASSOCIATED(list(ifield)%r3dptr) )THEN
                    ! Allow for varying size of third dimension
-                   r3dcount = r3dcount + SIZE(list(ifield)%r3dptr, index_z)
+                   r3dcount = r3dcount + nrecvp(irecv,1)
                 ELSE IF( ASSOCIATED(list(ifield)%i3dptr) )THEN
                    ! Allow for varying size of third dimension
-                   i3dcount = i3dcount + SIZE(list(ifield)%i3dptr, index_z)
+                   i3dcount = i3dcount + nrecvp(irecv,1)
                 END IF
              END IF
@@ -2957 +2992 @@
 
           IF ( r2dcount > 0 .OR. r3dcount > 0 ) THEN
-             CALL MPI_irecv (recvBuff(1,irecv),((r2dcount+r3dcount)*nrecvp(irecv,1)),     &
+             CALL MPI_irecv (recvBuff(1,irecv),(r2dcount+r3dcount),     &
                              MPI_DOUBLE_PRECISION, source(irecv), tag, mpi_comm_opa, &
                              exch_flags(handle,irecv,indexr), ierr)
           END IF
           IF ( i2dcount > 0 .OR. i3dcount > 0 ) THEN
-             CALL MPI_irecv (recvIBuff(1,irecv),((i2dcount+i3dcount)*nrecvp(irecv,1)),       &
+             CALL MPI_irecv (recvIBuff(1,irecv),(i2dcount+i3dcount),       &
                              MPI_INTEGER, source(irecv),tag, mpi_comm_opa, &
                              exch_flags(handle,irecv,indexr),ierr)
@@ -2993 +3028 @@
 
     ierr = 0
-    newSize = jpkdta*maxsendpts*nfields
+    newSize = maxsendpts*nfields
     IF( have_real_field .AND. newSize > sendBuffSize)THEN
        sendBuffSize=newSize
@@ -3010 +3045 @@
 
     IF (ierr .ne. 0) THEN
-       WRITE(*,*) 'ARPDBG: failed to allocate send buf'
        CALL ctl_stop('STOP','exchs_generic_list: unable to allocate send buff')
     END IF
@@ -3111 +3145 @@
                       DO j=jbeg, jend, 1
                          DO i=ibeg, iend, 1
-                            DO k=1, SIZE(list(ifield)%r3dptr, index_z), 1
-#else
-                      DO k=1, SIZE(list(ifield)%r3dptr, index_z), 1
+                            !DO k=1, SIZE(list(ifield)%r3dptr, index_z), 1
+                            DO k=1, nzsendp(ipatch,isend,1), 1
+#else
+                      !DO k=1, SIZE(list(ifield)%r3dptr, index_z), 1
+                      DO k=1, nzsendp(ipatch,isend,1), 1
                          DO j=jbeg, jend, 1
                             DO i=ibeg, iend, 1
@@ -3124 +3160 @@
                    
                       npacked =  nxsendp(ipatch,isend,1) * &
-                                 nysendp(ipatch,isend,1)
-                      rstart   = rstart + npacked*SIZE(list(ifield)%r3dptr, index_z)
-                      r3dcount = r3dcount + npacked*SIZE(list(ifield)%r3dptr, index_z)
+                                 nysendp(ipatch,isend,1) * &
+                                 nzsendp(ipatch,isend,1)
+                      rstart   = rstart + npacked
+                      r3dcount = r3dcount + npacked
+
                    END DO pack_patches3r
 
@@ -3143 +3181 @@
                       DO j=jbeg, jend, 1
                          DO i=ibeg, iend, 1
-                            DO k=1, SIZE(list(ifield)%i3dptr, index_z),1
-#else
-                      DO k=1, SIZE(list(ifield)%i3dptr, index_z),1
+                            !DO k=1, SIZE(list(ifield)%i3dptr, index_z),1
+                            DO k=1, nzsendp(ipatch,isend,1), 1
+#else
+                      !DO k=1, SIZE(list(ifield)%i3dptr, index_z),1
+                      DO k=1, nzsendp(ipatch,isend,1), 1
                          DO j=jbeg, jend, 1
                             DO i=ibeg, iend, 1
@@ -3155 +3195 @@
                       END DO
 
-                      istart   = istart +  nxs*nys*SIZE(list(ifield)%i3dptr, index_z)
-                      i3dcount = i3dcount + nxs*nys*SIZE(list(ifield)%i3dptr, index_z)
+                      npacked = nxsendp(ipatch,isend,1)* &
+                                nysendp(ipatch,isend,1)* &
+                                nzsendp(ipatch,isend,1)
+                      istart   = istart +  npacked
+                      i3dcount = i3dcount + npacked
                    END DO pack_patches3i
 
@@ -3172 +3215 @@
           ! Now do the send(s) for all fields
           IF(r2dcount > 0 .OR. r3dcount > 0 )THEN
-             CALL MPI_Isend(sendBuff(1,isend),(r2dcount+r3dcount),MPI_DOUBLE_PRECISION, &
-                            destination(isend),tag,mpi_comm_opa, &
+             CALL MPI_Isend(sendBuff(1,isend),(r2dcount+r3dcount), &
+                            MPI_DOUBLE_PRECISION,                  &
+                            destination(isend),tag,mpi_comm_opa,   &
                             exch_flags(handle,isend,indexs),ierr)
           END IF
@@ -3228 +3272 @@
 
                 ! Increment starting index for next field data in buffer
-                rstart = rstart + nrecvp(irecv,1)
+                rstart = ic + 1 !rstart + nrecvp(irecv,1)
 
              ELSE IF ( ASSOCIATED(list(ifield)%i2dptr) ) THEN
@@ -3250 +3294 @@
 
                 ! Increment starting index for next field data in buffer
-                istart = istart + nrecvp(irecv,1)
+                istart = ic + 1 !istart + nrecvp(irecv,1)
 
              ELSE IF (ASSOCIATED(list(ifield)%r3dptr) ) THEN
@@ -3264 +3308 @@
                    DO j=jbeg, jend, 1
                       DO i=ibeg, iend, 1
-                         DO k=1,SIZE(list(ifield)%r3dptr, index_z), 1
-#else
-                   DO k=1,  SIZE(list(ifield)%r3dptr, index_z), 1
+                         DO k=1, nzrecvp(ipatch,irecv,1), 1
+#else
+                   DO k=1, nzrecvp(ipatch,irecv,1), 1
                       DO j=jbeg, jend, 1
                          DO i=ibeg, iend, 1
@@ -3278 +3322 @@
 
                 ! Increment starting index for next field data in buffer
-                rstart = rstart + nrecvp(irecv,1)*SIZE(list(ifield)%r3dptr,index_z)
+                rstart = ic + 1 ! rstart + nrecvp(irecv,1) !*SIZE(list(ifield)%r3dptr,index_z)
 
              ELSE IF ( ASSOCIATED(list(ifield)%i3dptr) ) THEN
@@ -3292 +3336 @@
                    DO j=jbeg, jend, 1
                       DO i=ibeg, iend, 1
-                         DO k=1,SIZE(list(ifield)%i3dptr,index_z),1
-#else
-                   DO k=1,SIZE(list(ifield)%i3dptr,index_z),1
+                         DO k=1,nzrecvp(ipatch,irecv,1),1
+#else
+                   DO k=1,nzrecvp(ipatch,irecv,1),1
                       DO j=jbeg, jend, 1
                          DO i=ibeg, iend, 1
@@ -3306 +3350 @@
 
                 ! Increment starting index for next field data in buffer
-                istart = istart + nrecvp(irecv,1)*SIZE(list(ifield)%i3dptr,index_z)
+                istart = ic + 1 !istart + nrecvp(irecv,1) !*SIZE(list(ifield)%i3dptr,index_z)
 
              END IF
@@ -3395 +3439 @@
     CALL free_exch_handle(handle)
 
-    CALL prof_region_end(ARPEXCHS_LIST, iprofStat)
+    !CALL prof_region_end(ARPEXCHS_LIST, iprofStat)
     
   END SUBROUTINE exchs_generic_list
@@ -3431 +3475 @@
     ! *******************************************************************
     USE par_oce,     ONLY: wp, jpreci, jprecj, jpni, jpkdta
-    USE mapcomm_mod, ONLY: Iplus, Jplus, Iminus, Jminus, IplusJplus, &
-                           IminusJminus, IplusJminus, IminusJplus,   &
-                           nrecv, nsend, nrecvp, nsendp, nxsend,nysend,&
-                           destination,dirsend, dirrecv, &
-                           isrcsend, jsrcsend, idesrecv, jdesrecv, &
-                           isrcsendp,jsrcsendp,idesrecvp,jdesrecvp, &
-                           nxrecv,nyrecv,source, iesub, jesub, &
-                           MaxCommDir, MaxComm, idessend, jdessend, &
-                           nxsendp, nysendp, nxrecvp, nyrecvp,      &
-                           npatchsend, npatchrecv, &
-                           cyclic_bc
-    USE lib_mpp,     ONLY: mpi_comm_opa, ctl_stop
+    USE mapcomm_mod, ONLY: Iplus, Jplus, Iminus, Jminus, IplusJplus,   &
+                           IminusJminus, IplusJminus, IminusJplus,     &
+                           nrecv, nsend, nrecvp, nsendp,               &
+                           nrecvp2d, nsendp2d,  nxsend, nysend,        &
+                           destination,dirsend, dirrecv,               &
+                           isrcsend, jsrcsend, idesrecv, jdesrecv,     &
+                           isrcsendp,jsrcsendp,idesrecvp,jdesrecvp,    &
+                           nxrecv,source, iesub, jesub,         &
+                           MaxCommDir, MaxComm, idessend, jdessend,    &
+                           nxsendp, nysendp, nzsendp,                  &
+                           nxrecvp, nyrecvp, nzrecvp,                  &
+                           npatchsend, npatchrecv, cyclic_bc
+    USE lib_mpp,     ONLY: ctl_stop
+#if defined key_mpp_mpi
+    USE lib_mpp,     ONLY: mpi_comm_opa
+#endif
     USE dom_oce,     ONLY: narea
     USE in_out_manager, ONLY: numout
@@ -3474 +3522 @@
     INTEGER :: index  ! To hold index returned from MPI_waitany
     INTEGER, DIMENSION(3) :: isubsizes, istarts ! isizes
+#if defined key_mpp_mpi
     INTEGER :: status(MPI_status_size)
     INTEGER :: astatus(MPI_status_size,MaxComm)
+#endif
     LOGICAL, SAVE :: first_time = .TRUE.
 #if defined key_z_first
@@ -3489 +3539 @@
 
     !CALL prof_region_begin(ARPEXCHS_GENERIC, "Exchs_indiv", iprofStat)
-!    CALL timing_start('exchs_generic')
+    !CALL timing_start('exchs_generic')
 
     ierr = 0
 
-    ! Find out the sizes of the arrays.
-
-    kdim1 = 1
-    IF ( PRESENT(b3) ) THEN
-       kdim1 = SIZE(b3,dim=index_z)
-    ELSEIF ( PRESENT(ib3) ) THEN
-       kdim1 = SIZE(ib3,dim=index_z)
-    ELSEIF ( PRESENT(b2) ) THEN
-       kdim1 = SIZE(b2,dim=2)
-    ELSEIF ( PRESENT(ib2) ) THEN
-       kdim1 = SIZE(ib2,dim=2)
-    ENDIF
-
     ! Check nhexch is in range.
 
     IF ( nhexch.GT.jpreci ) THEN
-       STOP 'exchs: halo width greater than maximum'
+       CALL ctl_stop('STOP','exchs: halo width greater than maximum')
     ENDIF
 
@@ -3544 +3581 @@
        IF(.NOT. ALLOCATED(sendBuff))THEN
           ! Only allocate the sendBuff once
-          ALLOCATE(recvBuff(jpkdta*maxrecvpts,nrecv), &
-                   sendBuff(jpkdta*maxsendpts,nsend),stat=ierr)
+          ALLOCATE(recvBuff(maxrecvpts,nrecv), &
+                   sendBuff(maxsendpts,nsend),stat=ierr)
        ELSE
-          ALLOCATE(recvBuff(jpkdta*maxrecvpts,nrecv),stat=ierr)
+          ALLOCATE(recvBuff(maxrecvpts,nrecv),stat=ierr)
        END IF
     ELSE IF(PRESENT(ib2) .OR. PRESENT(ib3))THEN
        IF(.NOT. ALLOCATED(sendIBuff))THEN
-          ALLOCATE(recvIBuff(jpkdta*maxrecvpts,nrecv), &
-                   sendIBuff(jpkdta*maxsendpts,nsend),stat=ierr)
+          ALLOCATE(recvIBuff(maxrecvpts,nrecv), &
+                   sendIBuff(maxsendpts,nsend),stat=ierr)
        ELSE
-          ALLOCATE(recvIBuff(jpkdta*maxrecvpts,nrecv),stat=ierr)
+          ALLOCATE(recvIBuff(maxrecvpts,nrecv),stat=ierr)
        END IF
     END IF
@@ -3578 +3615 @@
           !          that isn't used
           IF ( PRESENT(b2) ) THEN
-             CALL MPI_irecv (recvBuff(1,irecv),nrecvp(irecv,1),   &
+             CALL MPI_irecv (recvBuff(1,irecv),nrecvp2d(irecv,1), &
                              MPI_DOUBLE_PRECISION, source(irecv), &
                              tag, mpi_comm_opa,                   &
                              exch_flags(handle,irecv,indexr), ierr)
           ELSEIF ( PRESENT(ib2) ) THEN
+             CALL MPI_irecv (recvIBuff(1,irecv),nrecvp2d(irecv,1), &
+                             MPI_INTEGER, source(irecv),         &
+                             tag, mpi_comm_opa,                  &
+                             exch_flags(handle,irecv,indexr),ierr)
+          ELSEIF ( PRESENT(b3) ) THEN
+             CALL MPI_irecv (recvBuff(1,irecv),nrecvp(irecv,1),   &
+                             MPI_DOUBLE_PRECISION, source(irecv), &
+                             tag, mpi_comm_opa,                   &
+                             exch_flags(handle,irecv,indexr),ierr)
+          ELSEIF ( PRESENT(ib3) ) THEN
              CALL MPI_irecv (recvIBuff(1,irecv),nrecvp(irecv,1), &
                              MPI_INTEGER, source(irecv),         &
                              tag, mpi_comm_opa,                  &
                              exch_flags(handle,irecv,indexr),ierr)
-          ELSEIF ( PRESENT(b3) ) THEN
-             CALL MPI_irecv (recvBuff(1,irecv),nrecvp(irecv,1)*kdim1,   &
-                             MPI_DOUBLE_PRECISION, source(irecv), &
-                             tag, mpi_comm_opa,                   &
-                             exch_flags(handle,irecv,indexr),ierr)
-          ELSEIF ( PRESENT(ib3) ) THEN
-             CALL MPI_irecv (recvIBuff(1,irecv),nrecvp(irecv,1)*kdim1, &
-                             MPI_INTEGER, source(irecv),         &
-                             tag, mpi_comm_opa,                  &
-                             exch_flags(handle,irecv,indexr),ierr)
           ENDIF
-          IF ( ierr.NE.0 ) THEN
-             WRITE (numout,*) 'ARPDBG - irecv hit error'
-             CALL flush(numout)
-             CALL MPI_abort(mpi_comm_opa,1,ierr)
-          END IF
+          ! No point checking for MPI errors because default MPI error handler
+          ! aborts run without returning control to calling program.
+          !IF ( ierr.NE.0 ) THEN
+          !   WRITE (numout,*) 'ARPDBG - irecv hit error'
+          !   CALL flush(numout)
+          !   CALL MPI_abort(mpi_comm_opa,1,ierr)
+          !END IF
 
 #if defined DEBUG_COMMS
           WRITE (*,FMT="(I4,': exchs post recv : hand = ',I2,' dirn = ',I1,' src = ',I3,' tag = ',I4,' npoints = ',I6)") &
                   narea-1,handle,dirrecv(irecv), &
-                  source(irecv), tag, nrecvp(irecv,1)*kdim1
+                  source(irecv), tag, nrecvp(irecv,1)
 #endif
 
@@ -3635 +3674 @@
 
        IF ( enabled(dirsend(isend)) .AND. &
-            destination(isend).GE.0 .AND. nxsend(isend).GT.0 ) THEN
+            destination(isend) >= 0 .AND. nxsend(isend) > 0 ) THEN
 
           isrc = isrcsend(isend)
@@ -3647 +3686 @@
           IF(PRESENT(b3))THEN
              WRITE (*,FMT="(I4,': handle ',I4,' tag ',I4,' sending to ',I4,' data ',I4,' direction ',I3)") &  
-               narea-1, handle, tag, destination(isend),nsendp(isend,1)*kdim1,dirsend(isend)
+               narea-1, handle, tag, destination(isend),nsendp(isend,1),dirsend(isend)
           ELSE IF(PRESENT(b2))THEN
              WRITE (*,FMT="(I4,': handle ',I4,' tag ',I4,' sending to ',I4,' data ',I4,' direction ',I3)") &  
-               narea-1, handle, tag, destination(isend),nsendp(isend,1),dirsend(isend)
+               narea-1, handle, tag, destination(isend),nsendp2d(isend,1),dirsend(isend)
           END IF
 #endif
@@ -3673 +3712 @@
                 END DO
 
+!!$                ! For 'stupid' compiler that refuses to do a memcpy for above
 !!$                CALL do_real8_copy( nxsendp(patch,isend,1)*nysendp(patch,isend,1), &
 !!$                                    b2(istart,jstart),                             &
@@ -3709 +3749 @@
           ELSEIF ( PRESENT(b3) )THEN
 
-!             CALL timing_start('3dr_pack')
+             ! CALL timing_start('3dr_pack')
              ic = 0
              pack_patches3r: DO ipatch=1,npatchsend(isend,1)
@@ -3720 +3760 @@
                 DO j=jstart, jend, 1
                    DO i=istart, iend, 1
-                      DO k=1,kdim1,1
-#else
-                DO k=1,kdim1,1
+                      DO k=1,nzsendp(ipatch,isend,1),1
+#else
+                DO k=1,nzsendp(ipatch,isend,1),1
                    DO j=jstart, jend, 1
                       DO i=istart, iend, 1
@@ -3732 +3772 @@
                 END DO
              END DO pack_patches3r
-!             CALL timing_stop('3dr_pack')
+
+             ! CALL timing_stop('3dr_pack')
 
              CALL MPI_Isend(sendBuff(1,isend),ic,                  &
@@ -3740 +3781 @@
 
 #if defined DEBUG_COMMS
-          WRITE (*,FMT="(I4,': Isend of ',I3,' patches, ',I6,' points, to ',I3)") &
+             WRITE (*,FMT="(I4,': Isend of ',I3,' patches, ',I6,' points, to ',I3)") &
                      narea-1, npatchsend(isend,1),ic, &
                      destination(isend)
@@ -3756 +3797 @@
                  DO j=jstart, jend, 1
                     DO i=istart, iend, 1
-                       DO k=1,kdim1,1
-#else
-                 DO k=1,kdim1,1
+                       DO k=1,nzsendp(ipatch,isend,1),1
+#else
+                 DO k=1,nzsendp(ipatch,isend,1),1
                     DO j=jstart, jend, 1
                        DO i=istart, iend, 1
@@ -3775 +3816 @@
           ENDIF
 
-          IF ( ierr.NE.0 ) CALL MPI_abort(mpi_comm_opa,1,ierr)
+          !IF ( ierr.NE.0 ) CALL MPI_abort(mpi_comm_opa,1,ierr)
 
        ELSE
@@ -3785 +3826 @@
     ENDDO ! Loop over sends
 
-!    CALL timing_stop('mpi_sends')
+    ! CALL timing_stop('mpi_sends')
 
 #if ( defined DEBUG && defined DEBUG_EXCHANGE ) || defined DEBUG_COMMS
@@ -3793 +3834 @@
     ! Wait on the receives that were posted earlier
 
-!    CALL timing_start('mpi_recvs')
+    ! CALL timing_start('mpi_recvs')
 
     ! Copy just the set of flags we're interested in for passing 
@@ -3814 +3855 @@
           WRITE (*,"(I3,': ERROR: exchs_generic: MPI_waitany returned unrecognised error')") narea-1
        END IF
-       CALL ctl_stop('STOP')
+       CALL ctl_stop('STOP','exchs_generic: MPI_waitany returned error')
     END IF
 
@@ -3821 +3862 @@
           IF ( PRESENT(b2) ) THEN
 
-!             CALL timing_start('2dr_unpack')
+             ! CALL timing_start('2dr_unpack')
 
              ! Copy received data back into array
@@ -3839 +3880 @@
              END DO unpack_patches2r
 
-!             CALL timing_stop('2dr_unpack')
+             ! CALL timing_stop('2dr_unpack')
 
           ELSE IF ( PRESENT(ib2) ) THEN
@@ -3861 +3902 @@
            ELSE IF (PRESENT(b3) ) THEN
 
-!             CALL timing_start('3dr_unpack')
+              ! CALL timing_start('3dr_unpack')
              ic = 0
              unpack_patches3r: DO ipatch=1,npatchrecv(irecv,nhexch)
@@ -3872 +3913 @@
                 DO j=jstart, jend, 1
                    DO i=istart, iend, 1
-                      DO k=1,kdim1,1
-#else
-                DO k=1,kdim1,1
+                      DO k=1,nzrecvp(ipatch,irecv,1),1
+#else
+                DO k=1,nzrecvp(ipatch,irecv,1),1
                    DO j=jstart, jend, 1
                       DO i=istart, iend, 1
@@ -3881 +3922 @@
                          b3(i,j,k) = recvBuff(ic,irecv)
                       END DO
+#if defined key_z_first
+                      ! ARPDBG - wipe anything below the ocean bottom
+                      DO k=nzrecvp(ipatch,irecv,1)+1,jpk,1
+                         b3(i,j,k) = 0.0_wp
+                      END DO
+#endif
                    END DO
                 END DO
+
+                ! ARPDBG - wipe anything below the ocean bottom
+#if ! defined key_z_first
+                DO k=nzrecvp(ipatch,irecv,1)+1,jpk,1
+                   DO j=jstart, jend, 1
+                      DO i=istart, iend, 1
+                         b3(i,j,k) = 0.0_wp
+                      END DO
+                   END DO
+                END DO
+#endif 
+
              END DO unpack_patches3r
 
@@ -3899 +3958 @@
                 DO j=jstart, jend, 1
                    DO i=istart, iend, 1
-                      DO k=1,kdim1,1
-#else
-                DO k=1,kdim1,1
+                      DO k=1,nzrecvp(ipatch,irecv,1),1
+#else
+                DO k=1,nzrecvp(ipatch,irecv,1),1
                    DO j=jstart, jend, 1
                       DO i=istart, iend, 1
@@ -3915 +3974 @@
 
        CALL MPI_waitany (nrecv, exch_flags1d, irecv, status, ierr)
-       IF ( ierr.NE.0 ) CALL MPI_abort(mpi_comm_opa,1,ierr)
+       !IF ( ierr.NE.0 ) CALL MPI_abort(mpi_comm_opa,1,ierr)
 
     END DO ! while irecv != MPI_UNDEFINED
 
-!    CALL timing_stop('mpi_recvs')
+    ! CALL timing_stop('mpi_recvs')
 
     ! All receives done and unpacked so can deallocate the associated
     ! buffers
-    IF(ALLOCATED(recvBuff ))DEALLOCATE(recvBuff)
-    IF(ALLOCATED(recvIBuff))DEALLOCATE(recvIBuff)
+    !IF(ALLOCATED(recvBuff ))DEALLOCATE(recvBuff)
+    !IF(ALLOCATED(recvIBuff))DEALLOCATE(recvIBuff)
 
 #if defined DEBUG_COMMS
@@ -3940 +3999 @@
     !          loop!
     IF ( cyclic_bc .AND. (jpni.EQ.1) ) THEN
+
+       ! Find out the sizes of the arrays.
+       kdim1 = 1
+       IF ( PRESENT(b3) ) THEN
+          kdim1 = SIZE(b3,dim=index_z)
+       ELSEIF ( PRESENT(ib3) ) THEN
+          kdim1 = SIZE(ib3,dim=index_z)
+       ENDIF
+
 
        IF ( enabled(Iplus) ) THEN
@@ -3996 +4064 @@
        ENDIF
 
-    ENDIF
+    ENDIF ! cyclic_bc .AND. jpni == 1
 
     ! Copy just the set of flags we're interested in for passing to  
@@ -4009 +4077 @@
     IF( ALLOCATED(recvIBuff) )DEALLOCATE(recvIBuff)
 
-!    CALL timing_stop('exchs_generic')
+    ! CALL timing_stop('exchs_generic')
     !CALL prof_region_end(ARPEXCHS_GENERIC, iprofStat)
 
@@ -4313 +4381 @@
     CALL prof_region_begin(ARPNORTHLISTCOMMS, "NorthList", iprofStat)
 
+#if defined key_mpp_mpi
+
     ! If we get into this routine it's because : North fold condition and mpp 
     ! with more than one PE across i : we deal only with the North condition
 
     ! Set no. of rows from a module parameter that is also used in exchtestmod
+    ! and mpp_ini_north
     ijpj = num_nfold_rows
 
@@ -5244 +5315 @@
     CALL prof_region_end(NORTHLISTSCATTER, iprofStat)
 
+#endif /* key_mpp_mpi */
+
     CALL prof_region_end(ARPNORTHLISTCOMMS, iprofStat)
 
@@ -5255 +5328 @@
     !!
     !! ** Purpose :
-    !!      Ensure proper north fold horizontal bondary condition in mpp configuration
-    !!      in case of jpn1 > 1 (for 2d array )
+    !!      Ensure proper north fold horizontal bondary condition in mpp 
+    !!      configuration in case of jpn1 > 1 (for 2d array )
     !!
     !! ** Method :
@@ -5266 +5339 @@
     !!   8.5  !  03-09  (J.M. Molines ) For mpp folding condition at north
     !!                                  from lbc routine
-    !!   9.0  !  03-12  (J.M. Molines ) encapsulation into lib_mpp, coding rules of lbc_lnk
+    !!   9.0  !  03-12  (J.M. Molines ) encapsulation into lib_mpp, coding 
+    !!                                  rules of lbc_lnk
     !!----------------------------------------------------------------------
     USE par_oce,     ONLY : jpni, jpi, jpj
@@ -5273 +5347 @@
     USE mapcomm_mod, ONLY : pielb, piesub
     USE lib_mpp,     ONLY : ctl_stop
+    USE arpdebugging, ONLY: dump_array
     IMPLICIT none
     !! * Arguments
@@ -5287 +5362 @@
     !! * Local declarations
 
-    INTEGER, PARAMETER :: ijpj = 4
+    INTEGER :: ijpj
     INTEGER :: ji, jj,  jr, jproc
     INTEGER :: ierr
@@ -5303 +5378 @@
     ! with more than one PE across i : we deal only with the North condition
 
+    ! Set local from public PARAMETER
+    ijpj = num_nfold_rows
+
     CALL prof_region_begin(ARPNORTHCOMMS2D, "North2D", iprofStat)
 
+#if defined key_mpp_mpi
+
     IF(.not. ALLOCATED(ztab2))THEN
 
-       ALLOCATE(ztab2(jpiglo,4),                &
-                znorthgloio2(nwidthmax,4,jpni), &
-                znorthloc2(nwidthmax,4),        &
+       ALLOCATE(ztab2(jpiglo,ijpj),                &
+                znorthgloio2(nwidthmax,ijpj,ndim_rank_north), &
+                znorthloc2(nwidthmax,ijpj),        &
                 STAT=ierr)
        IF(ierr .ne. 0)THEN
@@ -5321 +5401 @@
     ijpjm1=ijpj-1
 
-    ! put the last 4 jlines of pt2d into znorthloc2
+    ! put the last ijpj jlines of pt2d into znorthloc2
     znorthloc2(:,:) = 0_wp ! because of padding for nwidthmax
-    DO ij = 1, ijpj, 1
+
+    ! jeub is the upper j limit of current domain in global coords
+    !
+    !                      |======================= jpjglo     ^ 
+    !    <Trimmed>         |                                  /|\
+    !                      |----------------------- jpjglo-1   |
+    !                      |                                   |
+    ! |---------jeub--------------------------------           
+    ! |                    |                                   j
+    ! |--------------------------------------------            
+    ! |                    |                                   |
+    ! |--------------------------------------------            |
+    !
+    ! No. of trimmed rows = jpjglo - jeub
+    ! No. of valid rows for n-fold = ijpj - <no. trimmed rows>
+    !                              = ijpj - jpjglo + jeub
+    ! Need an iterator that ends with max value ijpj and has (ijpj-jpjglo+jeub)
+    ! distinct values so start point must be:
+    !  ij_start = ijpj - (ijpj-jpjglo+jeub) + 1 = jpjglo - jeub + 1
+    ! => if jeub == jpjglo then we recover a starting value of 1.
+    !    if jeub == jpjglo - 10 then ij_start = 11 so no loop iterations
+    !    will be performed.
+
+#if defined NO_NFOLD_GATHER
+    ! Post receives for other PE's north-fold data
+    DO iproc = 1, ndim_rank_north, 1
+
+       IF( iproc-1 ==  nrank_north(iproc) ) CYCLE ! Skip this PE
+
+       CALL MPI_IRecv(znorthgloio2(), north_pts(iproc), MPI_DOUBLE_PRECISION, &
+                      nrank_north(iproc), iproc, tag, ncomm_north,            &
+                      nexch_flag(iproc) )
+    END DO
+#endif
+
+    DO ij = jpjglo - jeub + 1, ijpj, 1
+
        jj = nlcj - ijpj + ij
        znorthloc2(nldi:nlei,ij)=pt2d(nldi:nlei,jj)
     END DO
+
+!    CALL dump_array(0,'znorthloc2',znorthloc2,withHalos=.TRUE.,toGlobal=.FALSE.)
 
     IF (npolj /= 0 ) THEN
@@ -5335 +5453 @@
                        znorthgloio2,itaille,MPI_DOUBLE_PRECISION,  &
                        0, ncomm_north, ierr)
+
     ENDIF
 
     IF (narea == north_root+1 ) THEN
        ! recover the global north array
+       ! ztab2 has full width of global domain
        ztab2(:,:) = 0_wp
 
@@ -5350 +5470 @@
        END DO
 
+!       CALL dump_array(0,'ztab2',ztab2,withHalos=.TRUE.,toGlobal=.FALSE.)
 
        ! 2. North-Fold boundary conditions
@@ -5495 +5616 @@
       ENDIF
 
-      ! put in the last ijpj jlines of pt2d znorthloc2
-      DO ij = 1, ijpj, 1
+      ! Put the last ijpj jlines of pt2d into znorthloc2 while allowing
+      ! for any trimming of domain (see earlier comments and diagram)
+      DO ij = jpjglo - jeub + 1, ijpj, 1
          jj = nlcj - ijpj + ij
          pt2d(nldi:nlei,jj)= znorthloc2(nldi:nlei,ij)
       END DO
+
+#endif /* key_mpp_mpi */
 
       CALL prof_region_end(ARPNORTHCOMMS2D, iprofStat)
@@ -5512 +5636 @@
     !!
     !! ** Purpose :
-    !!      Ensure proper north fold horizontal bondary condition in mpp configuration
-    !!      in case of jpn1 > 1 (for 2d array )
+    !!      Ensure proper north fold horizontal bondary condition in mpp 
+    !!      configuration in case of jpn1 > 1 (for 2d array )
     !!
     !! ** Method :
@@ -5545 +5669 @@
     !! * Local declarations
 
-    INTEGER, PARAMETER :: ijpj = 4
+    INTEGER :: ijpj
     INTEGER :: ji, jj,  jr, jproc
     INTEGER :: ierr
@@ -5561 +5685 @@
     ! with more than one PE across i : we deal only with the North condition
 
+#if defined key_mpp_mpi
+
+    ijpj = num_nfold_rows
+    ijpjm1=ijpj - 1
+
+
      IF(.not. ALLOCATED(ztab2))THEN
 
-        ALLOCATE(ztab2(jpiglo,4),                &
-                 znorthgloio2(nwidthmax,4,jpni), &
-                 znorthloc2(nwidthmax,4),        &
+        ALLOCATE(ztab2(jpiglo,ijpj),                &
+                 znorthgloio2(nwidthmax,ijpj,jpni), &
+                 znorthloc2(nwidthmax,ijpj),        &
                  STAT=ierr)
         IF(ierr .ne. 0)THEN
@@ -5575 +5705 @@
     ! ---------------
 
-    ijpjm1=ijpj - 1
-
-    ! put in znorthloc2 the last 4 jlines of ib2
+    ! Put the last ijpj jlines of ib2 into znorthloc2 while allowing
+    ! for any trimming of domain (see earlier comments and diagram in
+    ! mpp_lbc_north_2d).
     znorthloc2(:,:) = 0  ! because of padding for nwidthmax
-    DO ij = 1, ijpj, 1
+    DO ij = jpjglo - jeub + 1, ijpj, 1
        jj = nlcj - ijpj + ij
        znorthloc2(nldi:nlei,ij)=ib2(nldi:nlei,jj)
@@ -5602 +5732 @@
           ilei=nleit (jproc)
           iilb=pielb(jproc)
-          WRITE (*,*)'ARPDBG, jproc = ',jproc,' ildi, ilei, iilb and ijpj = ',ildi, ilei, iilb, ijpj
+          !WRITE (*,*)'ARPDBG, jproc = ',jproc,' ildi, ilei, iilb and ijpj = ',&
+          !            ildi, ilei, iilb, ijpj
           ztab2(iilb:iilb+piesub(jproc)-1,1:ijpj) = &
                                      znorthgloio2(ildi:ilei,1:ijpj,jr)
@@ -5740 +5871 @@
             ilei=nleit (jproc)
             iilb=pielb(jproc)
-            znorthgloio2(ildi:ilei,1:ijpj,jr)=ztab2(iilb:iilb+piesub(jproc)-1,1:ijpj)
+            znorthgloio2(ildi:ilei,1:ijpj,jr) = &
+                                ztab2(iilb:iilb+piesub(jproc)-1,1:ijpj)
          END DO
 
@@ -5752 +5884 @@
       ENDIF
 
-      ! put in the last ijpj jlines of ib2 znorthloc2
-      DO ij = 1, ijpj, 1
+      ! put in the last ijpj jlines of ib2 from znorthloc2 while allowing
+      ! for any trimming of domain (see earlier comments and diagram in
+      ! mpp_lbc_north_2d).
+      DO ij = jpjglo - jeub + 1, ijpj, 1
          jj = nlcj - ijpj + ij
          ib2(nldi:nlei,jj)= znorthloc2(nldi:nlei,ij)
       END DO
       WRITE(*,*) 'ARPDBG: finished in mpp_lbc_north_i2d'
+
+#endif /* key_mpp_mpi */
+
    END SUBROUTINE mpp_lbc_north_i2d
 
@@ -5797 +5934 @@
 
      !! * Local declarations
-     INTEGER, PARAMETER :: ijpj = 4
+     INTEGER :: ijpj
      INTEGER :: ji, jj, jk, jr, jproc
      INTEGER :: ierr
@@ -5814 +5951 @@
      ! mpp with more than one proc across i : we deal only with the North 
      ! condition
+#if defined key_mpp_mpi
+
+     ijpj = num_nfold_rows
+     ijpjm1=ijpj - 1
 
      IF(.not. ALLOCATED(ztab))THEN
 
-        ALLOCATE(ztab(jpiglo,4,jpk),                &
-                 znorthgloio(nwidthmax,4,jpk,jpni), &
-                 znorthloc(nwidthmax,4,jpk),        &
+        ALLOCATE(ztab(jpiglo,ijpj,jpk),                &
+                 znorthgloio(nwidthmax,ijpj,jpk,jpni), &
+                 znorthloc(nwidthmax,ijpj,jpk),        &
                  STAT=ierr)
         IF(ierr .ne. 0)THEN
@@ -5835 +5976 @@
      ! ---------------
 
-    ijpjm1=ijpj - 1
-
-    ! Put the last ijpj jlines of pt3d into znorthloc
-    !ARPDBG znorthloc(:,:,:) = 999_wp ! because of padding for nwidthmax - 999 is
-                              ! for debugging
-#if defined key_z_first
-    DO ij = 1, ijpj, 1
+    ! Put the last ijpj jlines of pt3d into znorthloc while allowing
+    ! for any trimming of domain (see earlier comments and diagram in
+    ! mpp_lbc_north_2d).
+    ! Have to initialise all to zero in case chunks are missing due to domain
+    ! trimming
+    znorthloc(:,:,:) = 0.0_wp
+#if defined key_z_first
+    DO ij = jpjglo - jeub + 1, ijpj, 1
        jj = nlcj - ijpj + ij
        DO jk = 1, jpk 
 #else
     DO jk = 1, jpk 
-       DO ij = 1, ijpj, 1
+       DO ij = jpjglo - jeub + 1, ijpj, 1
           jj = nlcj - ijpj + ij
 #endif
@@ -5856 +5998 @@
     IF (npolj /= 0 ) THEN
        ! Build in proc 0 of ncomm_north the znorthgloio
-       !ARPDBG znorthgloio(:,:,:,:) = 0_wp
 
 #ifdef key_mpp_shmem
@@ -5875 +6016 @@
     IF (narea == north_root+1 ) THEN
        ! recover the global north array
-       !ARPDBG ztab(:,:,:) = 0_wp
+       ztab(:,:,:) = 0_wp
 
        DO jr = 1, ndim_rank_north
@@ -5890 +6031 @@
        ! ===============
 #if defined key_z_first
-
 
        ! 2. North-Fold boundary conditions
@@ -6179 +6319 @@
 #endif
 
-    ! put in the last ijpj jlines of pt3d znorthloc
-#if defined key_z_first
-    DO ij = 1, ijpj, 1
+    ! put in the last ijpj jlines of pt3d znorthloc while allowing
+    ! for any trimming of domain (see earlier comments and diagram in
+    ! mpp_lbc_north_2d).
+#if defined key_z_first
+    DO ij = jpjglo - jeub + 1, ijpj, 1
        jj = nlcj - ijpj + ij
        DO jk = 1 , jpk 
 #else
     DO jk = 1 , jpk 
-       DO ij = 1, ijpj, 1
+       DO ij = jpjglo - jeub + 1, ijpj, 1
           jj = nlcj - ijpj + ij
 #endif
@@ -6194 +6336 @@
 
     CALL prof_region_end(NORTH3DSCATTER, iprofStat)
+
+#endif /* key_mpp_mpi */
 
   END SUBROUTINE mpp_lbc_north_3d
@@ -6235 +6379 @@
 
      !! * Local declarations
-     INTEGER, PARAMETER :: ijpj = 4
-     INTEGER, PARAMETER :: ijpjm1 = ijpj - 1
+     INTEGER :: ijpj
+     INTEGER :: ijpjm1
      INTEGER :: ii, ji, jj, jk, jr, jproc
      INTEGER :: ierr
@@ -6254 +6398 @@
      ! mpp with more than one proc across i : we deal only with the North 
      ! condition
+
+     ijpj = num_nfold_rows
+     ijpjm1 = ijpj - 1
 
      IF(.not. ALLOCATED(ztab))THEN
@@ -6269 +6416 @@
      ! ---------------
 
-    ! put in znorthloc the last ijpj jlines of pt3d
-    znorthloc(:,:,:) = 0 ! because of padding for nwidthmax
-#if defined key_z_first
-    DO ij = 1, ijpj, 1
+    ! put in znorthloc the last ijpj jlines of pt3d while allowing
+    ! for any trimming of domain (see earlier comments and diagram in
+    ! mpp_lbc_north_2d).
+    znorthloc(:,:,:) = 0 ! because of padding for nwidthmax and domain
+                         ! trimming
+#if defined key_z_first
+    DO ij = jpjglo - jeub + 1, ijpj, 1
        jj = nlcj - ijpj + ij
        DO jk = 1, jpk 
 #else
     DO jk = 1, jpk 
-       DO ij = 1, ijpj, 1
+       DO ij = jpjglo - jeub + 1, ijpj, 1
           jj = nlcj - ijpj + ij
 #endif
@@ -6608 +6758 @@
 #endif
 
-    ! put in the last ijpj jlines of pt3d znorthloc
-#if defined key_z_first
-    DO ij = 1, ijpj, 1
+    ! put in the last ijpj jlines of pt3d znorthloc while allowing
+    ! for any trimming of domain (see earlier comments and diagram in
+    ! mpp_lbc_north_2d).
+#if defined key_z_first
+    DO ij = jpjglo - jeub + 1, ijpj, 1
        jj = nlcj - ijpj + ij
        DO ii = nldi, nlei, 1
@@ -6616 +6768 @@
 #else
     DO jk = 1 , jpk 
-       DO ij = 1, ijpj, 1
+       DO ij = jpjglo - jeub + 1, ijpj, 1
           jj = nlcj - ijpj + ij
           DO ii = nldi, nlei, 1

branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/LBC/exchtestmod.F90

@@ -7 +7 @@
                                                          ! once these tests are complete
 #if defined key_mpp_rkpart
-                        compare_exch_methods= .TRUE., & ! 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       = .TRUE.   ! Whether to use the halo
+                        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     = .FALSE., &
+                        use_exch_list       = .FALSE.   ! Whether to use the halo
                                                         ! packing API for the tests
                                                         ! - NOT currently working!
 #else
-                        compare_exch_methods= .FALSE., & ! Do both packed and 
-                        do_integer_tests    = .FALSE., & ! indiv. exchanges and compare output
+                        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     = .TRUE.,  &
+                        do_real_timings     = .FALSE., &
                         use_exch_list       = .FALSE.   ! Whether to use the halo
                                                         ! packing API for the tests
@@ -42 +43 @@
   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
@@ -62 +66 @@
   !====================================================================
 
-  SUBROUTINE mpp_test_comms(depth)
+  SUBROUTINE mpp_test_comms(depth, lmaxdepth)
     USE par_oce,      ONLY: jpi, jpj, jpk, jpreci
     USE par_kind,     ONLY: wp
     USE mapcomm_mod,  ONLY: Iminus, Iplus, Jminus, Jplus, NONE
-    USE lib_mpp,      ONLY: mpi_comm_opa, ctl_stop
+    USE lib_mpp,      ONLY: ctl_stop, ctl_warn
     USE dom_oce,      ONLY: narea
-    USE arpdebugging, Only: dump_array
     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(:,:) :: depth ! Mask (1 for ocean, 0 for land)
+    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()
 
@@ -87 +99 @@
     ALLOCATE(r3d(jpi,jpj,jpk), r3d_2(jpi,jpj,jpk), r3d_3(jpi,jpj,jpk), &
              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
@@ -135 +151 @@
        IF(narea == 1) WRITE (*,*) 'Test 3 done.'
 
-       ! 4. Test halo exchanges for a 3D REAL array at 'T' point...
+       ! 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'
@@ -144 +169 @@
        CALL exch_test(name, 'T', GLOBAL_LOCN_TEST, depth, r3d1=r3d)
 
-       IF(narea == 1) WRITE (*,*) 'Test 4 done.'
-
-       ! 5. Test halo exchanges for a 3D REAL array at 'V' point...
+       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'
@@ -152 +177 @@
             isgn=-999, lfill=.TRUE.)
 
-       IF(narea == 1) WRITE (*,*) 'Test 5 done.'
-
-       ! 6. Test halo exchanges for _two_ 2D REAL arrays at 'V' point...
+       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'
@@ -160 +185 @@
             isgn=-999, lfill=.TRUE.)
 
-       IF(narea == 1) WRITE (*,*) 'Test 6 done.'
-
-       ! 7. Test halo exchanges for _two_ 3D REAL arrays at 'V' point...
+       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'
@@ -168 +193 @@
             isgn=-999, lfill=.TRUE.)
 
-       IF(narea == 1) WRITE (*,*) 'Test 7 done.'
-
-       ! 8. Test halo exchanges for _three_ 3D REAL arrays at 'T' point...
+       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'
@@ -177 +202 @@
                       isgn=-999, lfill=.TRUE.)
 
-       IF(narea == 1) WRITE (*,*) 'Test 8 done.'
+       IF(narea == 1) WRITE (*,*) 'Test 9 done.'
 
     END IF
@@ -302 +327 @@
     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
@@ -316 +343 @@
        ! Generate a timing report
        CALL timing_finalize()
-       ! Dirty way of killing NEMO
+       ! Dirty way of causing NEMO to stop immediately
        CALL ctl_stop('STOP', 'Stopping now that comms tests are complete')
     END IF
@@ -331 +358 @@
     USE exchmod,      ONLY: add_exch, bound_exch_list, bound_exch
     USE lbclnk,       ONLY: lbc_lnk
-    USE arpdebugging, ONLY: dump_array
+    USE lib_mpp,      ONLY: ctl_warn
     USE dom_oce,      ONLY: narea
     IMPLICIT none
@@ -360 +387 @@
     !!-----------------------------------------------------------------------
 
+#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
@@ -1356 +1388 @@
                          lfill, stat)
     USE par_kind,     ONLY: wp
-    USE par_oce,      ONLY: jpi, jpj, jpk, jpreci
+    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
@@ -1374 +1406 @@
     LOGICAL,                          INTENT(out)          :: stat
     ! Locals
-    INTEGER :: ik, ij, ii
+    INTEGER :: ik, ij, ii, ipt
     LOGICAL :: hit_error, local_lfill
     INTEGER :: gVal, jstart, jstop, istart, istop
@@ -1495 +1527 @@
              DO ii=1,jpi,1
 #endif
-                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(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)
+                           narea-1, gridType, ik, ii, ij, r3d(ii, ij, ik)
 #else
-                          narea-1, gridType, ii, ij, ik, r3d(ii, ij, ik)
-#endif
+                           narea-1, gridType, ii, ij, ik, r3d(ii, ij, ik)
+#endif
+                   END IF
                 END IF
              END DO
@@ -1517 +1552 @@
                 DO ii=istart,istop,1
 
-                   gval = gcoords_to_int(ii,ij,ik)
-
-                   ! 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(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, &
@@ -1530 +1588 @@
 #endif
                            INT(r3d(ii, ij, ik)), gVal
-                      WRITE (LOG_UNIT,"(I4,': depth(',I3,',',I3,') = ',I2)") &
+                         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)
-                      hit_error = .TRUE.
+                            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
@@ -1556 +1618 @@
        END DO
 
-!!$          IF( ANY( MASK=(r2d(istart:istop,jstart:jstop) < 0.0) ) )THEN
-!!$             hit_error = .TRUE.
-!!$
-!!$             DO ij=jstart, jstop, 1
-!!$                DO ii=istart, istop, 1
-!!$                   IF(r2d(ii,ij) < 0.0)THEN
-!!$                      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
-!!$          END IF
-
        IF(testType .eq. GLOBAL_LOCN_TEST)THEN
 
           DO ij=jstart, jstop, 1
              DO ii=istart, istop, 1
-                gval = gcoords_to_int(ii,ij)
+
+                      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
@@ -1595 +1667 @@
              DO ij=jstart,jstop,1
                 DO ii=1,jpi,1
-                   IF(i3d(ii,ij,ik) < 0.0)THEN
+                   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
@@ -1618 +1692 @@
                 DO ii=istart,istop,1
                    gval = gcoords_to_int(ii,ij,ik)
-                   IF(i3d(ii,ij,ik) /= gval )THEN
+                   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
@@ -1677 +1752 @@
   END SUBROUTINE array_check
 
-  FUNCTION gcoords_to_int(ii, ij, ik) RESULT(value)
+  FUNCTION gcoords_to_int(ii, ij, ik, are_global) RESULT(value)
      USE dom_oce,      ONLY: nimpp, njmpp, nldi, nldj
      IMPLICIT None
-     ! Encode the specified global coordinates into a single
-     ! floating point number.
+     ! 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
-
-     value = (ii + nimpp - 1)*1000000 + &
-             (ij + njmpp - 1)*1000
+     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
@@ -1695 +1786 @@
   END FUNCTION gcoords_to_int
 
+
 END MODULE exchtestmod

branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/LBC/lbclnk.F90

@@ -18 +18 @@
    USE exchmod          ! Comms for irregular domain decomposition
 
+   ! This is important - it determines which set of comms routines are
+   ! called when lbc_lnk() is invoked.
    INTERFACE lbc_lnk
 #if  defined key_mpp_rkpart
@@ -104 +106 @@
 
 
-   SUBROUTINE lbc_lnk_3d( pt3d, cd_type, psgn, cd_mpp, pval )
+   SUBROUTINE lbc_lnk_3d( pt3d, cd_type, psgn, cd_mpp, pval, lzero )
+      USE lib_mpp, ONLY: ctl_stop
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE lbc_lnk_3d  ***
@@ -124 +127 @@
       CHARACTER(len=3)                , INTENT(in   ), OPTIONAL ::   cd_mpp    ! MPP only (here do nothing)
       REAL(wp)                        , INTENT(in   ), OPTIONAL ::   pval      ! background value (for closed boundaries)
+      LOGICAL                         , INTENT(in   ), OPTIONAL ::   lzero     ! Whether to zero halos on closed boundaries
+
       !!
       REAL(wp) ::   zland
@@ -132 +137 @@
       ENDIF
 
+      IF( PRESENT( lzero ) )THEN
+         CALL ctl_stop('STOP','lbc_lnk_3d: IMPLEMENT lzero option!')
+      ENDIF
 
       IF( PRESENT( cd_mpp ) ) THEN
@@ -197 +205 @@
 
 
-   SUBROUTINE lbc_lnk_2d( pt2d, cd_type, psgn, cd_mpp, pval )
+   SUBROUTINE lbc_lnk_2d( pt2d, cd_type, psgn, cd_mpp, pval, lzero )
+      USE lib_mpp, ONLY: ctl_stop
       !!---------------------------------------------------------------------
       !!                 ***  ROUTINE lbc_lnk_2d  ***
@@ -214 +223 @@
       CHARACTER(len=3)            , INTENT(in   ), OPTIONAL ::   cd_mpp    ! MPP only (here do nothing)
       REAL(wp)                    , INTENT(in   ), OPTIONAL ::   pval      ! background value (for closed boundaries)
+      LOGICAL                     , INTENT(in   ), OPTIONAL ::   lzero    ! Whether to zero halos on closed boundaries
       !!
       REAL(wp) ::   zland
@@ -220 +230 @@
       IF( PRESENT( pval ) ) THEN   ;   zland = pval      ! set land value (zero by default)
       ELSE                         ;   zland = 0.e0
+      ENDIF
+
+      IF( PRESENT( lzero ) )THEN
+         CALL ctl_stop('STOP','lbc_lnk_2d: IMPLEMENT lzero option!')
       ENDIF
 

branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/LBC/lib_mpp.F90

@@ -279 +279 @@
       WRITE(ldtxt(ii),*) '      size in bytes of exported buffer   nn_buffer   = ', nn_buffer     ;   ii = ii + 1
       WRITE(ldtxt(ii),*) '      whether to trim dry points         nn_pttrim   = ', nn_pttrim     ;   ii = ii + 1
-      WRITE(ldtxt(ii),*) '      number of cores per compute node   nn_cpn      = ', nn_cpnode     ;   ii = ii + 1
+      WRITE(ldtxt(ii),*) '      number of cores per compute node   nn_cpnode   = ', nn_cpnode     ;   ii = ii + 1
 #if defined key_agrif
       IF( .NOT. Agrif_Root() ) THEN
@@ -469 +469 @@
 
 #if defined key_mpp_rkpart
-      CALL ctl_stop('mpp_lnk_3d: should not have been called when key_mpp_rkpart defined!')
+      CALL ctl_stop('STOP', &
+                    'mpp_lnk_3d: should not have been called when key_mpp_rkpart defined!')
       RETURN
 #endif
@@ -2308 +2309 @@
       !!----------------------------------------------------------------------
       !
+#if defined key_mpp_rkpart
+      WRITE(*,*)'ARPDBG - should not be calling this version of mpp_ini_north!'
+      CALL MPI_ABORT(mpi_comm_opa, -1)
+      RETURN
+#endif
+
       njmppmax = MAXVAL( njmppt )
       !

branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/LBC/mapcomm_mod.F90

@@ -1 +1 @@
 MODULE mapcomm_mod
   USE in_out_manager, ONLY: numout, lwp
-  USE par_oce,        ONLY: jpiglo, jpjglo, jpreci, jprecj, jpi
+  USE par_oce,        ONLY: jpiglo, jpjglo, jpreci, jprecj, jpi, jpk
   USE dom_oce,        ONLY: nlei, nlej, nldi, nldj, nlci, nlcj, narea, &
                             nleit, nlejt, nldit, nldjt, nlcit, nlcjt, &
@@ -7 +7 @@
   IMPLICIT NONE
 
-!#define ARPDEBUG
+#define ARPDEBUG
 
   PRIVATE
@@ -62 +62 @@
                                        isrcrecv,jsrcrecv, &
                                        idessend,jdessend, &
-                                       nxsend,nysend,     &
+                                       nxsend,nysend,nzsend, &
                                        idesrecv,jdesrecv, &
-                                       nxrecv,nyrecv
+                                       nxrecv,nyrecv,nzrecv
   INTEGER, SAVE :: nsend,nrecv
 
   ! SMP 22 Sep 2009
-  ! Alternate run-length encoded communications lists
+  ! Alternative, run-length encoded communications lists
   ! omitting permanently dry points.
-  ! Of these, idessendp, jdessendp, isrcrecp, jsrcrecvp
+  ! Of these, isrcrecp, jsrcrecvp
   ! are set up but not currently used,
   ! and could be eliminated.
@@ -79 +79 @@
   INTEGER, SAVE, DIMENSION(MaxPatch,MaxComm,jpreci) :: &
                                     isrcsendp, jsrcsendp,&
-                                    !idessendp, jdessendp,&
-                                    nxsendp, nysendp,    &
+                                    nxsendp, nysendp, nzsendp, &
                                     isrcrecvp, jsrcrecvp,&
                                     idesrecvp, jdesrecvp,&
-                                    nxrecvp, nyrecvp
+                                    nxrecvp, nyrecvp, nzrecvp
   INTEGER, SAVE, DIMENSION(MaxComm,jpreci) :: npatchsend, npatchrecv
   ! Total number of points in each message
-  INTEGER, SAVE, DIMENSION(MaxComm,jpreci) :: nsendp, nrecvp
+  INTEGER, SAVE, DIMENSION(MaxComm,jpreci) :: nsendp, nsendp2d, nrecvp, nrecvp2d
 
   ! Process dependent partitioning information.
@@ -132 +131 @@
                ,south = (/ 0, 0, 1, 0, 1, 0, 0, 1 /) &
                ,north = (/ 0, 0, 0, 1, 0, 1, 1, 0 /)
-!                          1  2  3  4  5  6  7  8
-!                          W  E  S  N SW NE NW  SE
+                         ! 1  2  3  4  5  6  7  8
+                         ! W  E  S  N SW NE NW  SE
 
   ! cyclic_bc     True if a cyclic boundary condition is to be applied
@@ -150 +149 @@
   INTEGER, PARAMETER :: LAND = 0
 
-  ! nextra is a safety factor because NEMO actually computes
-  ! its wet/dry mask in dommsk _after_ it has smoothed the
-  ! bathymetry read from file (when ln_sco is set). This means
-  ! that points on the coast that are dry here can actually 
-  ! subsequently become wet. Therefore, rather than trim to a point
-  ! immediately next to a wet point, we back off nextra points.
-  INTEGER, PARAMETER :: nextra = 2
+  ! Rather than trim to a point immediately next to a wet point, we 
+  ! back off nextra points. If we don't do this then the sea-ice
+  ! computation goes wrong because it does use values over the land
+  ! that immediately border the ocean.
+  INTEGER, SAVE :: nextra
 
   ! Public routines
@@ -164 +161 @@
   PUBLIC :: MaxComm,nsend,nrecv,nxsend,nysend,destination,dirrecv, &
             dirsend,isrcsend,jsrcsend,idesrecv, jdesrecv,          &
-            nxrecv,nyrecv,source, cyclic_bc, idessend, jdessend
-
-  PUBLIC :: nsendp,nrecvp,npatchsend,npatchrecv, &
-            nxsendp,nysendp, nxrecvp,nyrecvp,    &
+            nxrecv, nyrecv, source, cyclic_bc, idessend, jdessend
+
+  PUBLIC :: nsendp,nsendp2d,nrecvp,nrecvp2d,npatchsend,npatchrecv, &
+            nxsendp,nysendp,nzsendp,nxrecvp,nyrecvp,nzrecvp,       &
             idesrecvp,jdesrecvp,isrcsendp,jsrcsendp
 
@@ -190 +187 @@
   PUBLIC :: trimmed, nidx, eidx, sidx, widx, nextra
 
-  ! Switch for outputting px mapping to file
-  !LOGICAL, PARAMETER :: outmap = .TRUE.
-
   ! Switch for trimming dry points from halo swaps
-  LOGICAL, PARAMETER :: msgtrim = .TRUE.
+  LOGICAL, PARAMETER :: msgtrim   = .TRUE.
+
+  ! Switch for trimming points below ocean floor from halo swaps
+  !LOGICAL, PARAMETER :: msgtrim_z = .TRUE. ! .FALSE.
+  LOGICAL, PUBLIC, SAVE      :: msgtrim_z
 
 CONTAINS
 
-  SUBROUTINE mapcomms ( depth, nx, ny, jperio, ierr )
+  SUBROUTINE mapcomms ( depth, ibotlvl, nx, ny, jperio, ierr )
     !!------------------------------------------------------------------
     ! Maps out the communications requirements for the partitioned
@@ -209 +207 @@
     ! Subroutine arguments.
     INTEGER, INTENT(in) :: nx, ny
-    INTEGER, INTENT(in) :: depth(nx,ny)! Global mask: 0 for land, 1 for ocean
-    INTEGER, INTENT(in) :: jperio      ! Periodicity of the mesh
+    INTEGER, INTENT(in) :: depth(nx,ny)  ! Global mask: 0 for land, 1 for ocean
+    INTEGER, INTENT(in) :: ibotlvl(nx,ny)! Last vert level above sea floor
+    INTEGER, INTENT(in) :: jperio        ! Periodicity of the mesh
     INTEGER, INTENT(out):: ierr
 
@@ -268 +267 @@
     nxsend = -999
     nysend = -999
+    nzsend = -999
     dirrecv = -999
     source = -999
@@ -274 +274 @@
     nxrecv = -999
     nyrecv = -999
+    nzrecv = -999
 
     ! For each of the eight communication directions on a 2d grid of
@@ -457 +458 @@
           CALL addsend (nsend,Iplus,procid(iproc), &
                         isrcs,jsrcs,idess,jdess,   &
-                        nxs,nys,depth,ierr)
+                        nxs,nys,depth,ibotlvl,ierr)
           IF ( ierr.NE.0 ) RETURN
 
@@ -464 +465 @@
           CALL addrecv (nrecv,Iminus,procid(iproc), &
                         isrcr,jsrcr,idesr,jdesr,    &
-                        nxr,nyr,depth,ierr)
+                        nxr,nyr,depth,ibotlvl,ierr)
           IF ( ierr.NE.0 ) RETURN
 #if defined ARPDEBUG
@@ -641 +642 @@
 !         of border.
 
-          CALL addsend (nsend,Iminus,procid(iproc), &
-                        isrcs,jsrcs,idess,jdess,nxs,nys,depth,ierr)
+          CALL addsend (nsend,Iminus,procid(iproc),     &
+                        isrcs,jsrcs,idess,jdess,nxs,nys,&
+                        depth,ibotlvl,ierr)
           IF ( ierr.NE.0 ) RETURN
 #if defined ARPDEBUG
@@ -651 +653 @@
 #endif
 
-          CALL addrecv (nrecv,Iplus,procid(iproc), &
-                        isrcr,jsrcr,idesr,jdesr,nxr,nyr,depth,ierr)
+          CALL addrecv (nrecv,Iplus,procid(iproc),       &
+                        isrcr,jsrcr,idesr,jdesr,nxr,nyr, &
+                        depth,ibotlvl,ierr)
           IF ( ierr.NE.0 ) RETURN
 
@@ -853 +856 @@
 !         of border.
 
-          CALL addsend (nsend,Jplus,procid(iproc)      &
-                       ,isrcs,jsrcs,idess,jdess,nxs,nys,depth,ierr)
+          CALL addsend (nsend,Jplus,procid(iproc),       &
+                        isrcs,jsrcs,idess,jdess,nxs,nys, &
+                        depth,ibotlvl,ierr)
           IF ( ierr.NE.0 ) RETURN
 
-          CALL addrecv (nrecv,Jminus,procid(iproc) &
-                 ,isrcr,jsrcr,idesr,jdesr,nxr,nyr,depth,ierr)
+          CALL addrecv (nrecv,Jminus,procid(iproc),      &
+                        isrcr,jsrcr,idesr,jdesr,nxr,nyr, &
+                        depth,ibotlvl,ierr)
           IF ( ierr.NE.0 ) RETURN
 
@@ -1029 +1034 @@
 !         of border.
 
-          CALL addsend (nsend,Jminus,procid(iproc) &
-                       ,isrcs,jsrcs,idess,jdess,nxs,nys,depth,ierr)
+          CALL addsend (nsend,Jminus,procid(iproc),      &
+                        isrcs,jsrcs,idess,jdess,nxs,nys, &
+                        depth,ibotlvl,ierr)
           IF ( ierr.NE.0 ) RETURN
 
-          CALL addrecv (nrecv,Jplus,procid(iproc) &
-                 ,isrcr,jsrcr,idesr,jdesr,nxr,nyr,depth,ierr)
+          CALL addrecv (nrecv,Jplus,procid(iproc),       &
+                        isrcr,jsrcr,idesr,jdesr,nxr,nyr, &
+                        depth,ibotlvl,ierr)
           IF ( ierr.NE.0 ) RETURN
 
@@ -1057 +1064 @@
 
       ! Loop over the four corner directions
+      ! i = 1  2  3  4  5  6  7 8
+      !     W  E  S  N SW NE NW SE
+
 
       DO i=5,8
@@ -1063 +1073 @@
 
         addcorner = .FALSE.
-
-        ! Look at the processors to the X and Y directions.
-
-!!$        iprocx = iprocmap(ielb-west(i)+east(i)*iesub &
-!!$                         ,jelb+north(i)*(jesub-1))
-
-!!$        iprocy = iprocmap(ielb+east(i)*(iesub-1)     &
-!!$                         ,jelb-south(i)+north(i)*jesub)
 
         ! i1 is to be x-coord just OUTSIDE our domain
@@ -1080 +1082 @@
            i1 = ielb
            i2 = ielb
-           IF(ilbext)THEN
+           IF(ilbext .AND. (.NOT. trimmed(widx,narea)) )THEN
               i2 = i2+west(i) ! If on W boundary with cyclic bc's, ielb _is_ the halo column
                               ! so add 1 to move inside domain
@@ -1086 +1088 @@
               i1 = i1-west(i)
            END IF
-           IF(iubext)THEN
+           IF(iubext .AND. (.NOT. trimmed(eidx,narea)) )THEN
               ! If upper bound is on domain boundary then iesub already
               ! includes the halo column
@@ -1101 +1103 @@
         END IF
 
-! For a NW corner:
-!               | 
-!       iproc   |  iprocy
-!       ________|______
-!               |
-!       iprocx  |  Me
-!               |
+        ! For a NW corner:
+        !               | 
+        !       iproc   |  iprocy
+        !       ________|______
+        !               |
+        !       iprocx  |  Me
+        !               |
 
         ! x coord just OUTSIDE our domain but y INSIDE
@@ -1141 +1143 @@
 
              ! Ensure we don't include halos from the global borders if we
-              ! have cyclic E/W boundaries.
+             ! have cyclic E/W boundaries.
              ielb_iproc = pielb(iproc)
              ieub_iproc = pieub(iproc)
@@ -1178 +1180 @@
             ! Allow for wrap-around if necessary
             IF(cyclic_bc)THEN
-               IF(ldiff0 < 1) ldiff0 = ldiff0 + (jpiglo - 2) !ARPDBG -2 for consistency with procmap
-               IF(ldiff1 < 1) ldiff1 = ldiff1 + (jpiglo - 2) !ARPDBG -2 for consistency with procmap
+               IF(ldiff0 < 1)THEN
+                  !ARPDBG -2 for consistency with procmap
+                  ldiff0 = ldiff0 + (jpiglo - 2)
+               END IF
+               IF(ldiff1 < 1)THEN
+                  !ARPDBG -2 for consistency with procmap
+                  ldiff1 = ldiff1 + (jpiglo - 2)
+               END IF
             END IF
             nxs  (ihalo) = ihalo -  east(i)*(ldiff0-1) &
                                  -  west(i)*(ldiff1-1)
+            ! Have no cyclic b.c.'s in N/S direction so probably don't need
+            ! the following checks on ldiff{0,1}
             ldiff0 = pjelb(iprocc) - jeub
             IF(ldiff0 < 1) ldiff0 = ldiff0 + jpjglo
@@ -1194 +1204 @@
             isrcs(ihalo) = east(i) *(iesub-nxs(ihalo)) + nldi
             jsrcs(ihalo) = north(i)*(jesub-nys(ihalo)) + nldj
-            IF(cyclic_bc)THEN
+            IF( cyclic_bc )THEN
                IF( ilbext )THEN
                   ! nldi is still within halo for domains on W edge of
@@ -1221 +1231 @@
             ! Source for a receive must be in an internal region of the REMOTE domain
             isrcr(ihalo) = west(i)*(piesub(iprocc)-nxs(ihalo)) + nldit(iprocc)
-            IF(cyclic_bc)THEN
+            IF( cyclic_bc )THEN
+
+               ! This _could_ be a corner exchange wrapped around by the cyclic
+               ! boundary conditions:
+               !
+               !  ||------|                      ||
+               !  ||      |           |          ||
+               !  ||a_____|__ _ _     |          ||
+               !  ||                  -----------||
+               !  ||                    |       a||
+               !  ||                    |________||
+
                IF(pilbext(iprocc))THEN
                   ! nldi is still within halo for domains on E edge of
                   ! global domain
-                  isrcr(ihalo) = isrcr(ihalo) + 1
+                  isrcr(ihalo) = isrcr(ihalo) + east(i)
                ELSE IF(piubext(iprocc))THEN
                   ! Final column is actually halo for domains on W edge of
                   ! global domain
-                  isrcr(ihalo) = isrcr(ihalo) - 1
+                  isrcr(ihalo) = isrcr(ihalo) - west(i)
                END IF
             END IF
@@ -1281 +1302 @@
         IF ( addcorner ) THEN
 #if defined ARPDEBUG
-          WRITE (*,FMT="(I3,': ARPDBG adding corner send to ',I2,', dir = ',I1)") narea-1, procid(iprocc),i
-#endif
-          CALL addsend (nsend,i,procid(iprocc) &
-                       ,isrcs,jsrcs,idess,jdess,nxs,nys,depth,ierr)
+          WRITE (*,FMT="(I3,': ARPDBG adding corner send to ',I2,', dir = ',I1)") &
+                 narea-1, procid(iprocc),i
+#endif
+          CALL addsend (nsend,i,procid(iprocc),          &
+                        isrcs,jsrcs,idess,jdess,nxs,nys, &
+                        depth,ibotlvl,ierr)
           IF ( ierr.NE.0 ) RETURN
 
@@ -1291 +1314 @@
 
 #if defined ARPDEBUG
-          WRITE (*,FMT="(I3,': ARPDBG adding corner recv. from ',I3,', old dir = ',I1,' new dir = ',I1)") narea-1, procid(iprocc),i, j
-#endif
-          CALL addrecv (nrecv,j,procid(iprocc) &
-                       ,isrcr,jsrcr,idesr,jdesr,nxr,nyr,depth,ierr)
+          WRITE (*,FMT="(I3,': ARPDBG adding corner recv. from ',I3,', old dir = ',I1,' new dir = ',I1)") &
+                 narea-1, procid(iprocc),i, j
+#endif
+          CALL addrecv (nrecv,j,procid(iprocc),          &
+                        isrcr,jsrcr,idesr,jdesr,nxr,nyr, &
+                        depth,ibotlvl,ierr)
           IF ( ierr.NE.0 ) RETURN
 
@@ -1366 +1391 @@
 
       SUBROUTINE addsend ( icomm, dir, proc, isrc, jsrc, &
-                           ides, jdes, nx, ny, depth, ierr )
+                           ides, jdes, nx, ny, depth, ibotlvl, ierr )
 !!------------------------------------------------------------------
 !     Adds a send communication specified by the parameters dir through 
@@ -1372 +1397 @@
 !     icomm points to the last entry and is incremented and returned 
 !     if successful.
-
+!
 !     icomm                   int   in/out    Location in comms list.
 !     dir                     int   input     Direction.
@@ -1384 +1409 @@
 !     depth                         input     Global mask, 0 for land, 1 for wet
 !     ierr                    int   output    Error flag.
-
+!
 !               Mike Ashworth, CLRC Daresbury Laboratory, March 1999
 !               Stephen Pickles, STFC Daresbury Laboratory
@@ -1397 +1422 @@
          ! Global mask: 0 for land, 1 for ocean
          INTEGER, DIMENSION(:,:),    INTENT( in  ) :: depth
+         INTEGER, DIMENSION(:,:),    INTENT( in  ) :: ibotlvl
          INTEGER,                    INTENT( out ) :: ierr
          INTEGER, DIMENSION(jpreci), INTENT( in  ) :: isrc, jsrc, &
                                                       ides, jdes, nx, ny
          ! Values of corresponding input arguments after clipping
-         INTEGER, DIMENSION(jpreci) :: cisrc,cjsrc,cides,cjdes,cnx,cny
+         INTEGER, DIMENSION(jpreci) :: cisrc,cjsrc,cides,cjdes,cnx,cny,cnz
          ! Run-length encoded versions corresponding to above
-         INTEGER, dimension(MaxPatch,jpreci) :: risrc,rjsrc,rides,rjdes,rnx,rny
+         INTEGER, DIMENSION(MaxPatch,jpreci) :: risrc,rjsrc,rides,rjdes,rnx,rny,rnz
          ! Number of patches in run-length encoded message
          INTEGER, DIMENSION(jpreci) :: npatches
          INTEGER :: ihalo, ipatch
+         INTEGER :: nsendp_untrimmedz ! How many pts we'd be sending without
+                                      ! trimming in z direction
          ! Whether there is still a message after clipping
          LOGICAL :: something_left
@@ -1421 +1449 @@
          ! Can the message be clipped ?
 
-         CALL clip_msg(depth, isrc, jsrc, ides, jdes, nx, ny, &
-                       cisrc,cjsrc,cides,cjdes,cnx,cny,       &
-                       risrc,rjsrc,rides,rjdes,rnx,rny,       &
+         CALL clip_msg(depth, ibotlvl,                       &
+                       isrc, jsrc, ides, jdes, nx, ny,       &
+                       cisrc,cjsrc,cides,cjdes,cnx,cny,cnz,  &
+                       risrc,rjsrc,rides,rjdes,rnx,rny,rnz,  &
                        npatches, something_left)
                            
@@ -1451 +1480 @@
             nxsend(icomm)      = cnx(1)
             nysend(icomm)      = cny(1)
+            IF(msgtrim_z)THEN
+               nzsend(icomm)   = cnz(1)
+            ELSE
+               nzsend(icomm)   = jpk
+            END IF
+
+            ! Zero count of untrimmed pts to send
+            nsendp_untrimmedz = 0
 
             ! Also set up the comms lists encoded as the start points and
             ! lengths of the contiguous runs of wet points.
             DO ihalo=1,jpreci
-               nsendp(icomm,ihalo) = 0
+
+               nsendp2d(icomm,ihalo)   = 0
+               nsendp(icomm,ihalo)     = 0
                npatchsend(icomm,ihalo) = npatches(ihalo)
+
                DO ipatch=1,npatches(ihalo)
+
                   isrcsendp(ipatch,icomm,ihalo) = risrc(ipatch,ihalo)
                   jsrcsendp(ipatch,icomm,ihalo) = rjsrc(ipatch,ihalo)
-                  !idessendp(ipatch,icomm,ihalo) = rides(ipatch,ihalo)
-                  !jdessendp(ipatch,icomm,ihalo) = rjdes(ipatch,ihalo)
+
                   nxsendp(ipatch,icomm,ihalo)   = rnx(ipatch,ihalo)
                   nysendp(ipatch,icomm,ihalo)   = rny(ipatch,ihalo)
+                  IF(msgtrim_z)THEN
+                     nzsendp(ipatch,icomm,ihalo)= rnz(ipatch,ihalo)
+                  ELSE
+                     nzsendp(ipatch,icomm,ihalo) = jpk
+                  END IF
+
                   ! Sum the no. of points to be sent over all
-                  ! patches
-                  nsendp(icomm,ihalo) = nsendp(icomm,ihalo) &
-                           + rnx(ipatch,ihalo)*rny(ipatch,ihalo)
+                  ! patches for both 2D-array halos and 3D-array halos
+                  nsendp2d(icomm,ihalo) = nsendp2d(icomm,ihalo) +      &
+                                          nxsendp(ipatch,icomm,ihalo)* &
+                                          nysendp(ipatch,icomm,ihalo)
+                  nsendp(icomm,ihalo) = nsendp(icomm,ihalo) +          &
+                                          nxsendp(ipatch,icomm,ihalo)* &
+                                          nysendp(ipatch,icomm,ihalo)* &
+                                          nzsendp(ipatch,icomm,ihalo)
+                  IF(msgtrim_z)THEN
+                     nsendp_untrimmedz = nsendp_untrimmedz +           &
+                                          nxsendp(ipatch,icomm,ihalo)* &
+                                          nysendp(ipatch,icomm,ihalo)* &
+                                          jpk
+                  END IF
                END DO
             END DO
 
 #if defined ARPDEBUG
-            WRITE (*,FMT="(I3,': ARPDBG adding SEND:')") narea-1  
-            WRITE (*,FMT="(I3,': ARPDBG: icomm = ',I2)") narea-1,icomm
-            WRITE (*,FMT="(I3,': ARPDBG:   dir = ',I2)") narea-1,dirsend(icomm)
-            WRITE (*,FMT="(I3,': ARPDBG:  proc = ',I3)") narea-1,destination(icomm)
-            WRITE (*,FMT="(I3,': ARPDBG:  isrc = ',I3)") narea-1,isrcsend(icomm)
-            WRITE (*,FMT="(I3,': ARPDBG:  jsrc = ',I3)") narea-1,jsrcsend(icomm)
-            WRITE (*,FMT="(I3,': ARPDBG:  ides = ',I3)") narea-1,idessend(icomm)
-            WRITE (*,FMT="(I3,': ARPDBG:  jdes = ',I3)") narea-1,jdessend(icomm)
-            WRITE (*,FMT="(I3,': ARPDBG:    nx = ',I3)") narea-1,nxsend(icomm)
-            WRITE (*,FMT="(I3,': ARPDBG:    ny = ',I3)") narea-1,nysend(icomm)
-            WRITE (*,FMT="(I3,': ARPDBG:npatch = ',I3)") narea-1,npatches(1)
+            WRITE (*,FMT="(I4,': ARPDBG adding SEND:')") narea-1  
+            WRITE (*,FMT="(I4,': ARPDBG: icomm = ',I2)") narea-1,icomm
+            WRITE (*,FMT="(I4,': ARPDBG:   dir = ',I2)") narea-1,dirsend(icomm)
+            WRITE (*,FMT="(I4,': ARPDBG:  proc = ',I4)") narea-1,destination(icomm)
+            WRITE (*,FMT="(I4,': ARPDBG:  isrc = ',I4)") narea-1,isrcsend(icomm)
+            WRITE (*,FMT="(I4,': ARPDBG:  jsrc = ',I4)") narea-1,jsrcsend(icomm)
+            WRITE (*,FMT="(I4,': ARPDBG:  ides = ',I4)") narea-1,idessend(icomm)
+            WRITE (*,FMT="(I4,': ARPDBG:  jdes = ',I4)") narea-1,jdessend(icomm)
+            WRITE (*,FMT="(I4,': ARPDBG:    nx = ',I4)") narea-1,nxsend(icomm)
+            WRITE (*,FMT="(I4,': ARPDBG:    ny = ',I4)") narea-1,nysend(icomm)
+            WRITE (*,FMT="(I4,': ARPDBG:    nz = ',I4)") narea-1,nzsend(icomm)
+            WRITE (*,FMT="(I4,': ARPDBG:npatch = ',I3)") narea-1,npatches(1)
  
             DO ipatch=1,npatches(1)
-               WRITE (*,FMT="(I3,': ARPDBG:  patch ',I2,': isrc = ',I3)") &
+               WRITE (*,FMT="(I4,': ARPDBG:  patch ',I2,': isrc = ',I4)") &
                                   narea-1,ipatch,isrcsendp(ipatch,icomm,1)
-               WRITE (*,FMT="(I3,': ARPDBG:  patch ',I2,': jsrc = ',I3)") &
+               WRITE (*,FMT="(I4,': ARPDBG:  patch ',I2,': jsrc = ',I4)") &
                                   narea-1,ipatch,jsrcsendp(ipatch,icomm,1)
-               !WRITE (*,FMT="(I3,': ARPDBG:  patch ',I2,': ides = ',I3)") &
-               !                   narea-1,ipatch,idessendp(ipatch,icomm,1)
-               !WRITE (*,FMT="(I3,': ARPDBG:  patch ',I2,': jdes = ',I3)") &
-               !                   narea-1,ipatch,jdessendp(ipatch,icomm,1)
-               WRITE (*,FMT="(I3,': ARPDBG:  patch ',I2,':   nx = ',I3)") &
+               WRITE (*,FMT="(I4,': ARPDBG:  patch ',I2,':   nx = ',I4)") &
                                   narea-1,ipatch,nxsendp(ipatch,icomm,1)  
-               WRITE (*,FMT="(I3,': ARPDBG:  patch ',I2,':   ny = ',I3)") &
+               WRITE (*,FMT="(I4,': ARPDBG:  patch ',I2,':   ny = ',I4)") &
                                   narea-1,ipatch,nysendp(ipatch,icomm,1)  
+               WRITE (*,FMT="(I4,': ARPDBG:  patch ',I2,':   nz = ',I4)") &
+                                  narea-1,ipatch,nzsendp(ipatch,icomm,1)  
             END DO
 
-            WRITE (*,FMT="(I3,': ARPDBG:nsendp = ',I4)") narea-1,nsendp(icomm,1)
-            WRITE (*,FMT="(I3,': ARPDBG SEND ends')")    narea-1
+            WRITE (*,FMT="(I4,': ARPDBG:nsendp = ',I4)") narea-1,nsendp(icomm,1)
+            IF(msgtrim_z)THEN
+               WRITE (*,FMT="(I4,': ARPDBG:nsendp WITHOUT z trim = ',I4)") &
+                                                       narea-1,nsendp_untrimmedz
+            END IF
+            WRITE (*,FMT="(I4,': ARPDBG SEND ends')")    narea-1
 #endif
 
@@ -1508 +1568 @@
 
     SUBROUTINE addrecv ( icomm, dir, proc, isrc, jsrc, &
-                         ides, jdes, nx, ny, depth, ierr )
-!!------------------------------------------------------------------
-!     Adds a recv communication specified by the parameters dir through 
-!     to ny to the recv communication list at the next position. 
-!     icomm points to the last entry and is incremented and returned 
-!     if successful.
-
-!     icomm                   int   in/out    Location in comms list.
-!     dir                     int   input     Direction.
-!     proc                    int   input     Process id.
-!     isrc                    int   input     X coordinate of source data.
-!     jsrc                    int   input     Y coordinate of source data.
-!     ides                    int   input     X coordinate of destination data.
-!     jdes                    int   input     Y coordinate of destination data.
-!     nx                      int   input     Size in X of data to be sent.
-!     ny                      int   input     Size in Y of data to be sent.
-!     ierr                    int   output    Error flag.
-
-!            Mike Ashworth, CLRC Daresbury Laboratory, March 1999
-!!------------------------------------------------------------------
+                         ides, jdes, nx, ny, depth, ibotlvl, ierr )
+      !!------------------------------------------------------------------
+      !   Adds a recv communication specified by the parameters dir through 
+      !   to ny to the recv communication list at the next position. 
+      !   icomm points to the last entry and is incremented and returned 
+      !   if successful.
+      !
+      !   icomm                   int   in/out    Location in comms list.
+      !   dir                     int   input     Direction.
+      !   proc                    int   input     Process id.
+      !   isrc                    int   input     X coordinate of source data.
+      !   jsrc                    int   input     Y coordinate of source data.
+      !   ides                    int   input     X coordinate of dest. data.
+      !   jdes                    int   input     Y coordinate of dest. data.
+      !   nx                      int   input     Size in X of data to be sent.
+      !   ny                      int   input     Size in Y of data to be sent.
+      !   ierr                    int   output    Error flag.
+      !
+      !          Mike Ashworth, CLRC Daresbury Laboratory, March 1999
+      !!------------------------------------------------------------------
       IMPLICIT NONE
 
@@ -1535 +1595 @@
       INTEGER,                 INTENT(out)   :: ierr
       INTEGER, DIMENSION(:,:), INTENT( in  ) :: depth
+      INTEGER, DIMENSION(:,:), INTENT( in  ) :: ibotlvl
       INTEGER, DIMENSION(jpreci)             :: isrc, jsrc, ides, jdes, nx, ny
 
@@ -1540 +1601 @@
 
       ! Values of corresponding input arguments after clipping
-      INTEGER, DIMENSION(jpreci) :: cisrc,cjsrc,cides,cjdes,cnx,cny
+      INTEGER, DIMENSION(jpreci) :: cisrc,cjsrc,cides,cjdes,cnx,cny,cnz
       ! Run-length encoded versions corresponding to above
-      INTEGER, dimension(MaxPatch,jpreci) :: risrc,rjsrc,rides,rjdes,rnx,rny
+      INTEGER, dimension(MaxPatch,jpreci) :: risrc,rjsrc,rides,rjdes,rnx,rny,rnz
       ! Number of patches in run-length encoded message
       INTEGER, DIMENSION(jpreci) :: npatches
@@ -1561 +1622 @@
       ! Can the message be clipped ?
 
-      CALL clip_msg(depth, ides, jdes, isrc, jsrc, nx, ny, &
-                          cides,cjdes,cisrc,cjsrc,cnx,cny, &
-                          rides,rjdes,risrc,rjsrc,rnx,rny, &
+      CALL clip_msg(depth, ibotlvl,                  &
+                    ides, jdes, isrc, jsrc, nx, ny,  &
+                    cides,cjdes,cisrc,cjsrc,cnx,cny,cnz, &
+                    rides,rjdes,risrc,rjsrc,rnx,rny,rnz, &
                     npatches, something_left)                           
 
@@ -1589 +1651 @@
          idesrecv(icomm) = cides(1)
          jdesrecv(icomm) = cjdes(1)
+
          nxrecv(icomm)   = cnx(1)
          nyrecv(icomm)   = cny(1)
+         IF(msgtrim_z)THEN
+            nzrecv(icomm)   = cnz(1)
+         ELSE
+            nzrecv(icomm)   = jpk
+         END IF
 
          DO ihalo=1,jpreci
-            nrecvp(icomm,ihalo) = 0
+
+            nrecvp2d(icomm,ihalo) = 0
+            nrecvp(icomm,ihalo)   = 0
             npatchrecv(icomm,ihalo) = npatches(ihalo)
+
             DO ipatch=1,npatches(ihalo)
                isrcrecvp(ipatch,icomm,ihalo) = risrc(ipatch,ihalo)
@@ -1602 +1673 @@
                nxrecvp(ipatch,icomm,ihalo)   = rnx(ipatch,ihalo)
                nyrecvp(ipatch,icomm,ihalo)   = rny(ipatch,ihalo)
+               IF(msgtrim_z)THEN
+                  nzrecvp(ipatch,icomm,ihalo) = rnz(ipatch,ihalo)
+               ELSE
+                  nzrecvp(ipatch,icomm,ihalo) = jpk
+               END IF
+
                ! Sum the no. of points to be received over all
                ! patches
-               nrecvp(icomm,ihalo) = nrecvp(icomm,ihalo) + &
-                                     rnx(ipatch,ihalo)*rny(ipatch,ihalo)
+               nrecvp2d(icomm,ihalo) = nrecvp2d(icomm,ihalo) +           &
+                                            nxrecvp(ipatch,icomm,ihalo)* &
+                                            nyrecvp(ipatch,icomm,ihalo)
+                    
+               nrecvp(icomm,ihalo) = nrecvp(icomm,ihalo) +               &
+                                            nxrecvp(ipatch,icomm,ihalo)* &
+                                            nyrecvp(ipatch,icomm,ihalo)* &
+                                            nzrecvp(ipatch,icomm,ihalo)
             END DO
          END DO
@@ -1613 +1696 @@
          WRITE (*,FMT="(I3,': ARPDBG: icomm = ',I2)") narea-1,icomm
          WRITE (*,FMT="(I3,': ARPDBG:   dir = ',I2)") narea-1,dir
-         WRITE (*,FMT="(I3,': ARPDBG:  proc = ',I3)") narea-1,proc
-         WRITE (*,FMT="(I3,': ARPDBG:  isrc = ',I3)") narea-1,cisrc(1)
-         WRITE (*,FMT="(I3,': ARPDBG:  jsrc = ',I3)") narea-1,cjsrc(1)
-         WRITE (*,FMT="(I3,': ARPDBG:  ides = ',I3)") narea-1,cides(1)
-         WRITE (*,FMT="(I3,': ARPDBG:  jdes = ',I3)") narea-1,cjdes(1)
-         WRITE (*,FMT="(I3,': ARPDBG:    nx = ',I3)") narea-1,cnx(1)
-         WRITE (*,FMT="(I3,': ARPDBG:    ny = ',I3)") narea-1,cny(1)
+         WRITE (*,FMT="(I3,': ARPDBG:  proc = ',I4)") narea-1,proc
+         WRITE (*,FMT="(I3,': ARPDBG:  isrc = ',I4)") narea-1,isrcrecv(icomm)
+         WRITE (*,FMT="(I3,': ARPDBG:  jsrc = ',I4)") narea-1,jsrcrecv(icomm)
+         WRITE (*,FMT="(I3,': ARPDBG:  ides = ',I4)") narea-1,idesrecv(icomm)
+         WRITE (*,FMT="(I3,': ARPDBG:  jdes = ',I4)") narea-1,jdesrecv(icomm)
+         WRITE (*,FMT="(I3,': ARPDBG:    nx = ',I4)") narea-1,nxrecv(icomm)
+         WRITE (*,FMT="(I3,': ARPDBG:    ny = ',I4)") narea-1,nyrecv(icomm)
+         WRITE (*,FMT="(I3,': ARPDBG:    nz = ',I4)") narea-1,nzrecv(icomm)
          WRITE (*,FMT="(I3,': ARPDBG:npatch = ',I3)") narea-1,npatches(1)
          DO ipatch=1,npatches(1)
-            WRITE (*,FMT="(I3,': ARPDBG:  patch ',I2,': isrc = ',I3)") &
+            WRITE (*,FMT="(I3,': ARPDBG:  patch ',I2,': isrc = ',I4)") &
                                   narea-1,ipatch,isrcrecvp(ipatch,icomm,1)
-            WRITE (*,FMT="(I3,': ARPDBG:  patch ',I2,': jsrc = ',I3)") &
+            WRITE (*,FMT="(I3,': ARPDBG:  patch ',I2,': jsrc = ',I4)") &
                                   narea-1,ipatch,jsrcrecvp(ipatch,icomm,1)
-            WRITE (*,FMT="(I3,': ARPDBG:  patch ',I2,': ides = ',I3)") &
+            WRITE (*,FMT="(I3,': ARPDBG:  patch ',I2,': ides = ',I4)") &
                                   narea-1,ipatch,idesrecvp(ipatch,icomm,1)
-            WRITE (*,FMT="(I3,': ARPDBG:  patch ',I2,': jdes = ',I3)") &
+            WRITE (*,FMT="(I3,': ARPDBG:  patch ',I2,': jdes = ',I4)") &
                                   narea-1,ipatch,jdesrecvp(ipatch,icomm,1)
-            WRITE (*,FMT="(I3,': ARPDBG:  patch ',I2,':   nx = ',I3)") &
+            WRITE (*,FMT="(I3,': ARPDBG:  patch ',I2,':   nx = ',I4)") &
                                   narea-1,ipatch,nxrecvp(ipatch,icomm,1)  
-            WRITE (*,FMT="(I3,': ARPDBG:  patch ',I2,':   ny = ',I3)") &
+            WRITE (*,FMT="(I3,': ARPDBG:  patch ',I2,':   ny = ',I4)") &
                                   narea-1,ipatch,nyrecvp(ipatch,icomm,1)  
+            WRITE (*,FMT="(I3,': ARPDBG:  patch ',I2,':   nz = ',I4)") &
+                                  narea-1,ipatch,nzrecvp(ipatch,icomm,1)  
          END DO
          WRITE (*,FMT="(I3,': ARPDBG:nrecvp = ',I4)") narea-1,nrecvp(icomm,1)
@@ -1644 +1730 @@
 
 
-    SUBROUTINE clip_msg(depth, iloc, jloc, irem, jrem, nx, ny, &
-                               ciloc,cjloc,cirem,cjrem,cnx,cny,&
-                               riloc,rjloc,rirem,rjrem,rnx,rny,&
+    SUBROUTINE clip_msg(depth, ibotlvl,                      &
+                        iloc, jloc, irem, jrem, nx, ny,      &
+                        ciloc,cjloc,cirem,cjrem,cnx,cny,cnz, &
+                        riloc,rjloc,rirem,rjrem,rnx,rny,rnz, &
                         npatches, something_left)                           
-!!------------------------------------------------------------------
-!
-!     Clip any exterior rows or columns that are permanently dry 
-!     from the message.
-!
-!     depth          real    input   Land/sea mask - global coords
-!     iloc           int     input   local  X coordinate of data start
-!     jloc           int     input   local  Y coordinate of data start
-!     irem           int     input   remote X coordinate of data
-!     jrem           int     input   remote Y coordinate of data
-!     nx             int     input   Size in X of data to be sent
-!     ny             int     input   Size in Y of data to be sent
-!     ciloc          int     output  As iloc, after clipping
-!     cjloc          int     output  As jloc, after clipping
-!     cirem          int     output  As irem, after clipping
-!     cjrem          int     output  As jrem, after clipping
-!     cnx            int     output  As nx, after clipping
-!     cny            int     output  As ny, after clipping
-!
-!     The run-length encoded versions split a message into one
-!     or more patches, leaving out permanently dry rows/columns
-!
-!     riloc          int     output  As iloc, run-length encoded
-!     rjloc          int     output  As jloc, run-length encoded
-!     rirem          int     output  As irem, run-length encoded
-!     rjrem          int     output  As jrem, run-length encoded
-!     rnx            int     output  As nx, run-length encoded
-!     rny            int     output  As ny, run-length encoded
-!     npatches       int     output  Number of patches
-!
-!     something_left logical output
-!
-!     Stephen Pickles, STFC Daresbury Laboratory, August 2009
-!     - Written
-!     Stephen Pickles, STFC Daresbury Laboratory, September 2009
-!     - Added run-length encoding
-!!------------------------------------------------------------------
+      !!------------------------------------------------------------------
+      !
+      !     Clip any exterior rows or columns that are permanently dry 
+      !     from the message. Also remove any vertical levels that are
+      !     beneath the ocean floor.
+      !
+      !     depth          int     input   Land/sea mask - global coords
+      !     ibotlvl        int     input   Index of the last vertical level 
+      !                                    above sea floor
+      !     iloc           int     input   local  X coordinate of data start
+      !     jloc           int     input   local  Y coordinate of data start
+      !     irem           int     input   remote X coordinate of data
+      !     jrem           int     input   remote Y coordinate of data
+      !     nx             int     input   Size in X of data to be sent
+      !     ny             int     input   Size in Y of data to be sent
+      !     ciloc          int     output  As iloc, after clipping
+      !     cjloc          int     output  As jloc, after clipping
+      !     cirem          int     output  As irem, after clipping
+      !     cjrem          int     output  As jrem, after clipping
+      !     cnx            int     output  As nx, after clipping
+      !     cny            int     output  As ny, after clipping
+      !
+      !     The run-length encoded versions split a message into one
+      !     or more patches, leaving out permanently dry rows/columns
+      !
+      !     riloc          int     output  As iloc, run-length encoded
+      !     rjloc          int     output  As jloc, run-length encoded
+      !     rirem          int     output  As irem, run-length encoded
+      !     rjrem          int     output  As jrem, run-length encoded
+      !     rnx            int     output  As nx, run-length encoded
+      !     rny            int     output  As ny, run-length encoded
+      !     rnz            int     output  Max depth (level) of this patch
+      !     npatches       int     output  Number of patches
+      !
+      !     something_left logical output
+      !
+      !     Stephen Pickles, STFC Daresbury Laboratory, August 2009
+      !     - Written
+      !     Stephen Pickles, STFC Daresbury Laboratory, September 2009
+      !     - Added run-length encoding
+      !     Andrew Porter, STFC Daresbury Laboratory, January 2013
+      !     - Added trimming of levels below sea floor
+      !!------------------------------------------------------------------
       USE dom_oce,      ONLY: nimpp, njmpp
       IMPLICIT none
       ! Subroutine arguments.
       INTEGER, DIMENSION(:,:), INTENT(in) :: depth ! Global mask (0 dry, 1 wet)
+      INTEGER, DIMENSION(:,:), INTENT(in) :: ibotlvl ! Bottom level of ocean
       INTEGER, DIMENSION(jpreci) ::  iloc, jloc, irem, jrem, nx, ny
-      INTEGER, DIMENSION(jpreci) :: ciloc,cjloc,cirem,cjrem,cnx,cny
-      INTEGER, DIMENSION(MaxPatch,jpreci) :: riloc,rjloc,rirem,rjrem,rnx,rny
+      INTEGER, DIMENSION(jpreci) :: ciloc,cjloc,cirem,cjrem,cnx,cny,cnz
+      INTEGER, DIMENSION(MaxPatch,jpreci) :: riloc,rjloc,rirem,rjrem,rnx,rny,rnz
       INTEGER, DIMENSION(jpreci), INTENT(out) :: npatches
-      LOGICAL, INTENT(out)        :: something_left
+      LOGICAL,                    INTENT(out) :: something_left
       ! Local variables.
       INTEGER :: h, i, j, patch
       LOGICAL :: all_dry
 
-      ! i, j limits of the halo patch, in local co-ordinates
+      ! i, j, k limits of the halo patch, in local co-ordinates
       ! These are set from input arguments, then updated as we trim
       INTEGER :: ilo, ihi, jlo, jhi
@@ -1708 +1802 @@
       cnx(:)      = nx(:)
       cny(:)      = ny(:)
+      cnz(:)      = jpk
       riloc(1,:)  = iloc(:)
       rjloc(1,:)  = jloc(:)
@@ -1714 +1809 @@
       rnx(1,:)    = nx(:)
       rny(1,:)    = ny(:)
+      rnz(:,:)    = jpk
       npatches(:) = 1
       something_left = .TRUE.
@@ -1730 +1826 @@
 
         ! Can any points along the left (low i) edge be trimmed?
-        left_edge: DO i=ilo, ihi
+        left_edge: DO i=ilo, ihi - nextra
           DO j=jlo, jhi
              ! depth is global mask, i and j are local coords
-            IF (depth(i+nimpp-1,j+njmpp-1) .NE. land) EXIT left_edge
+             ! ARPDBG - not sure that nextra needed below?
+            !IF (depth(i+nimpp-1+nextra,j+njmpp-1) .NE. LAND) EXIT left_edge
+            IF (depth(i+nimpp-1,j+njmpp-1) .NE. LAND) EXIT left_edge
           END DO
           ciloc(h) = ciloc(h) + 1
@@ -1743 +1841 @@
           cnx(h) = 0
           cny(h) = 0
+          cnz(h) = 0
           ciloc(h) = iloc(h)
           npatches(h) = 0
@@ -1757 +1856 @@
 
         ! Can any points along the right (high i) edge be trimmed?
-        right_edge: DO i=ihi, ilo, -1
+        right_edge: DO i=ihi, ilo + nextra, -1
           DO j=jlo, jhi
+!            IF (depth(i+nimpp-1-nextra,j+njmpp-1) .ne. land) exit right_edge
             IF (depth(i+nimpp-1,j+njmpp-1) .ne. land) exit right_edge
           END DO
@@ -1808 +1908 @@
            make_patches_x: DO WHILE (patch .lt. MaxPatch)
 
-              add_sea_cols: DO WHILE (i .lt. ihi)
-                 all_dry = .true.
-                 DO j=jlo, jhi
-                    IF (depth(i+nimpp-1,j+njmpp-1) .NE. land) THEN
-                       all_dry = .FALSE.
-                    END IF
-                 END DO
-                 IF (all_dry) EXIT add_sea_cols
-                 i = i+1
-                 rnx(patch,h) = rnx(patch,h) + 1
-              END DO add_sea_cols
+              IF(i == ihi)THEN
+
+                 rnx(patch,h) = 1
+              ELSE
+
+                 add_sea_cols: DO WHILE (i .lt. ihi)
+                    ! Check this strip in y to see whether all points are dry
+                    !               !all_dry = .TRUE.
+                    !IF( ANY( depth(ilo+nimpp-1:ihi+nimpp-2,j+njmpp-1) .NE. LAND ) )all_dry = .FALSE.
+                    !IF( ALL( depth(i+nimpp-1,jlo+njmpp-1:jhi+njmpp-1) == LAND ) )EXIT add_sea_cols
+
+                    all_dry = .TRUE.
+                    DO j=jlo, jhi
+                       IF (depth(i+nimpp-1,j+njmpp-1) .NE. land) THEN
+                          all_dry = .FALSE.
+                       END IF
+                    END DO
+                    IF (all_dry) EXIT add_sea_cols
+
+                    rnx(patch,h) = rnx(patch,h) + 1
+                    i = i+1
+                 END DO add_sea_cols
+              END IF
 
               ! This patch is now finished.
+
+              ! Store max depth of ocean bottom in this patch. riloc holds the starting
+              ! point of current patch in local coords.
+              ! riloc(patch,h) + nimpp - 1 is same point in global coords
+              ! End of patch is then at <start> + <length> - 1
+              rnz(patch,h) = MAXVAL(ibotlvl(riloc(patch,h)+nimpp-1:              &
+                                            riloc(patch,h)+rnx(patch,h)+nimpp-2, &
+                                            jlo+njmpp-1:jhi+njmpp-1) )
+
               ! Skip land cols before starting the next patch.
 
@@ -1849 +1970 @@
            ! Finish the last patch
            rnx(npatches(h),h) = ihi - riloc(npatches(h),h) + 1
+           rnz(npatches(h),h) = MAXVAL(ibotlvl(riloc(npatches(h),h)+nimpp-1:                    &
+                                               riloc(npatches(h),h)+rnx(npatches(h),h)+nimpp-2, &
+                                               jlo+njmpp-1:jhi+njmpp-1) )
 
         ELSE
@@ -1856 +1980 @@
            rnx(1,h) = cnx(h)
 
-          make_patches_y: do while (patch .lt. MaxPatch)
-
-            add_sea_rows: do while (j .lt. jhi)
-              all_dry = .true.
-              do i=ilo, ihi
-                if (depth(i+nimpp-1,j+njmpp-1) .ne. land) then
-                  all_dry = .false.
-                end if
-              end do
-              if (all_dry) exit add_sea_rows
-              j = j+1
-              rny(patch,h) = rny(patch,h) + 1
-            end do add_sea_rows
-
-            ! This patch is now finished.
-            ! Skip land rows before starting the next patch.
-
-            skip_land_rows: do while (j .lt. jhi)
-              do i=ilo, ihi
-                if (depth(i+nimpp-1,j+njmpp-1) .ne. land) then
-                  exit skip_land_rows
-                end if
-              end do
-              j = j+1
-            end do skip_land_rows
+           make_patches_y: DO WHILE (patch .lt. MaxPatch)
+
+              add_sea_rows: DO WHILE (j .lt. jhi)
+
+!                 IF( ALL( depth(ilo+nimpp-1:ihi+nimpp-1,                   &
+!                                j+njmpp-1)               == LAND ) )EXIT add_sea_rows
+
+                 all_dry = .TRUE.
+                 DO i=ilo, ihi
+                    if (depth(i+nimpp-1,j+njmpp-1) .ne. land) then
+                       all_dry = .FALSE.
+                    end if
+                 END DO
+                 IF (all_dry) EXIT add_sea_rows
+
+                 rny(patch,h) = rny(patch,h) + 1
+                 j = j+1
+              END DO add_sea_rows
+
+              ! This patch is now finished.
+
+              ! Store max depth of ocean bottom in this patch
+              rnz(patch,h) = MAXVAL(ibotlvl(ilo+nimpp-1:ihi+nimpp-1, &
+                                            rjloc(patch,h)+njmpp-1:  &
+                                            rjloc(patch,h)+rny(patch,h)+njmpp-2) )
+
+              ! Skip land rows before starting the next patch.
+
+              skip_land_rows: DO WHILE (j .lt. jhi)
+                 DO i=ilo, ihi
+                    IF (depth(i+nimpp-1,j+njmpp-1) .NE. LAND) THEN
+                       EXIT skip_land_rows
+                    END IF
+                 END DO
+                 j = j+1
+              END DO skip_land_rows
               
-            ! No more wet points?
-            if (j .ge. jhi) exit make_patches_y
-
-            ! Start next patch
-            patch = patch + 1
-            npatches(h) = patch
-            riloc(patch,h)  = ilo
-            rjloc(patch,h)  = j
-            rirem(patch,h)  = cirem(h)
-            rjrem(patch,h)  = cjrem(h)+j-jlo
-            rnx(patch,h)    = cnx(h)
-            rny(patch,h)    = 0
-
-          end do make_patches_y
-
-          ! Finish the last patch
-          rny(npatches(h),h) = jhi - rjloc(npatches(h),h) + 1
-
-        end if
+              ! No more wet points?
+              IF (j .ge. jhi) EXIT make_patches_y
+
+              ! Start next patch
+              patch = patch + 1
+              npatches(h) = patch
+              riloc(patch,h)  = ilo
+              rjloc(patch,h)  = j
+              rirem(patch,h)  = cirem(h)
+              rjrem(patch,h)  = cjrem(h)+j-jlo
+              rnx(patch,h)    = cnx(h)
+              rny(patch,h)    = 0
+
+           END DO make_patches_y
+
+           ! Finish the last patch
+           rny(npatches(h),h) = jhi - rjloc(npatches(h),h) + 1
+           rnz(npatches(h),h) = MAXVAL(ibotlvl(ilo+nimpp-1:ihi+nimpp-1,      &
+                                               rjloc(npatches(h),h)+njmpp-1: &
+                                               rjloc(npatches(h),h)+rny(npatches(h),h)+njmpp-2) )
+
+        END IF
+
+        ! Max depth for whole message is the maximum of the maximum depth of each
+        ! patch.
+        cnz(h) = MAXVAL(rnz(:,h))
 
       END DO haloes

branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/LBC/mppini.F90

@@ -136 +136 @@
       INTEGER, DIMENSION(jpni,jpnj) ::   iimppt, ijmppt, ilcit, ilcjt   ! local workspace
       INTEGER, ALLOCATABLE, DIMENSION(:,:) :: imask ! Local fake global land mask
+      INTEGER, ALLOCATABLE, DIMENSION(:,:), TARGET :: maxdepth ! Local fake global max depth mask
       !!----------------------------------------------------------------------
 
@@ -427 +428 @@
 
       ! ARPDBG - test comms setup
-      ALLOCATE(imask(jpiglo,jpjglo))
+      ALLOCATE(imask(jpiglo,jpjglo),maxdepth(jpiglo,jpjglo))
       imask(:,:) = 1
-      CALL mpp_test_comms(imask)
-      DEALLOCATE(imask)
+      maxdepth(:,:) = jpk
+      CALL mpp_test_comms(imask, maxdepth)
+      DEALLOCATE(imask, maxdepth)
 
    END SUBROUTINE mpp_init

branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/LBC/partition_mod.F90

@@ -1 +1 @@
 MODULE partition_mod
-   USE par_oce, ONLY: jpni, jpnj, jpnij, jpi, jpj, jpim1, jpjm1, jpij, &
+   USE par_oce, ONLY: jpni, jpnj, jpi, jpj, jpim1, jpjm1, jpij, &
                       jpreci, jprecj, jpk, jpkm1, jperio, jpiglo, jpjglo
    USE dom_oce, ONLY: ln_zco, nbondi, nbondj, nidom, npolj, &
@@ -15 +15 @@
                       nwidthmax,   &  ! Width of widest northern domain
                       narea           ! ID of local area (= rank + 1)
-   USE lib_mpp,        ONLY: mppsize, mppsync, mpi_comm_opa, mpp_ini_north, &
-                             ctl_stop, MAX_FACTORS, xfactors, yfactors,     &
-                             nn_pttrim, nn_cpnode
+#if defined key_mpp_mpi
+   USE lib_mpp,        ONLY: mppsize, mppsync, mpi_comm_opa,                &
+                             MAX_FACTORS, xfactors, yfactors, nn_pttrim,    &
+                             nn_cpnode
+#endif
+   USE lib_mpp,        ONLY: ctl_stop, ctl_warn
    USE in_out_manager, ONLY: numout, lwp
    USE mapcomm_mod,    ONLY: ielb, ieub, mapcomms, pielb, pjelb, pieub, pjeub,&
                              iesub, jesub, jeub, ilbext, iubext, jubext,      &
                              jlbext, pnactive, piesub, pjesub, jelb, pilbext, &
-                             piubext, pjlbext, pjubext,                       &
+                             piubext, pjlbext, pjubext, nextra,              &
                              nprocp   ! No. of PEs to partition over
    USE iom,            ONLY: wp, jpdom_unknown, iom_open, iom_get, iom_close
@@ -31 +34 @@
                                                  ! (1 for ocean, 0 for land)
                                                  ! set in nemogcm.F90
+   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:,:), TARGET :: ibotlevel ! Holds the bottom level of the ocean at each grid point - used for trimming halos in z direction
 
    ! Parameters for the cost function used when evaluating different 
@@ -67 +71 @@
    REAL(wp),PARAMETER :: pv_awful = 1.0e20
 
-! #define PARTIT_DEBUG
-
-   PUBLIC imask, smooth_bathy
+#define PARTIT_DEBUG
+
+   PUBLIC imask, ibotlevel, smooth_global_bathy, global_bot_level, partition_mask_alloc
    PUBLIC mpp_init3, partition_rk, partition_mca_rk, write_partition_map
 
 CONTAINS
+
+   SUBROUTINE partition_mask_alloc(xsize, ysize, ierr)
+      !!------------------------------------------------------------------
+      !!                  ***  ROUTINE partition_mask_alloc  ***
+      !!
+      !! Called from nemogcm to allocate the masks that are members of 
+      !! this module
+      !!
+      !!------------------------------------------------------------------
+      INTEGER, INTENT(in) :: xsize, ysize
+      INTEGER, INTENT(out):: ierr
+
+      ALLOCATE(imask(xsize,ysize), ibotlevel(xsize,ysize), Stat=ierr)
+
+   END SUBROUTINE partition_mask_alloc
+
 
    SUBROUTINE mpp_init3()
@@ -240 +260 @@
 
       ! Map out the communications for the partitioned domain.
-      CALL mapcomms (imask, jpiglo, jpjglo, jperio, ierr)
+      CALL mapcomms (imask, ibotlevel, jpiglo, jpjglo, jperio, ierr)
       IF ( ierr.NE.0 ) THEN
         IF ( lwp ) WRITE(numout,*) 'Communications mapping failed : ',ierr
@@ -247 +267 @@
 
       ! Prepare mpp north fold
-      IF (jperio >= 3 .AND. jperio <= 6 .AND. jpni > 1 ) THEN
-         CALL mpp_ini_north
-      END IF
+#if defined key_mpp_mpi
+      ! This invokes the version of the routine contained in this module
+      ! and not the original in lib_mpp.F90
+      CALL mpp_ini_north()
+#endif
 
 ! From mppini_2.h90:
@@ -271 +293 @@
 
       ! ARPDBG - test comms setup
-      CALL mpp_test_comms(imask)
+      CALL mpp_test_comms(imask, ibotlevel)
 
       ! Free array holding mask used for partitioning
@@ -326 +348 @@
       ENDIF
 
-      CALL flio_dom_set ( jpnij, narea-1, idid, iglo, iloc, iabsf, iabsl, &
+#if defined key_mpp_mpi
+      CALL flio_dom_set ( mppsize, narea-1, idid, iglo, iloc, iabsf, iabsl, &
                           ihals, ihale, 'BOX', nidom)
+#endif
 
     END SUBROUTINE mpp_init_ioipsl  
@@ -347 +371 @@
       USE iom,     ONLY: jpiglo, jpjglo, wp
       USE par_oce, ONLY: jpni, jpnj
+#if defined key_mpp_mpi
       USE lib_mpp, ONLY: mppsize
+#endif
       IMPLICIT NONE
 
@@ -355 +381 @@
       INTEGER, INTENT(in)        :: mask(:,:)
       ! Local variables
+#if defined key_mpp_mpi
       INTEGER, DIMENSION(MAX_FACTORS) :: fx,fy
+#endif
       INTEGER                    :: f,gnactive &
             ,i,ifax,ifin,ifx,ify,ilb,iproc,ist,isub,isub_old &
@@ -365 +393 @@
       ! Clear the error flag.
       ierr = 0
+
+#if defined key_mpp_mpi
 
       ! IMPORTANT: Set the number of PEs to partition over (mapcomm_mod 
@@ -398 +428 @@
       CALL finish_partition()
 
+#endif
+
     END SUBROUTINE partition_rk
 
@@ -404 +436 @@
 #if defined key_mpp_mpi
        USE mpi
-#endif
-       USE lib_mpp, ONLY: mppsize, ctl_stop, mpi_comm_opa, &
+       USE lib_mpp, ONLY: mppsize, mpi_comm_opa, &
                           nxfactors, nyfactors, xfactors, yfactors
+#endif
+       USE lib_mpp, ONLY: ctl_stop
        USE dom_oce, ONLY: narea
        IMPLICIT NONE
@@ -423 +456 @@
        ! Local variables
        INTEGER :: ii
+#if defined key_mpp_mpi
        INTEGER, DIMENSION(MAX_FACTORS) :: fx, fy, factors
        INTEGER, DIMENSION(MAX_FACTORS) :: df, multiplicity
+#endif
        INTEGER :: nfx, nfy, nfactors, ndf, nperms
        INTEGER :: check_nprocx, check_nprocy, check_nprocp
@@ -437 +472 @@
        INTEGER  :: best_perm
        REAL(wp), DIMENSION(2,pv_num_scores) :: best, gbest, wrst, gwrst
+
+#if defined key_mpp_mpi
 
        ! NEMO only has narea public and not the actual PE rank so
@@ -641 +678 @@
       END IF
 
+      ! Set corresponding NEMO variables for PE grid, even though it is now
+      ! rather irregular
+      jpni = nprocx
+      jpnj = nprocy
+
       IF (lwp) THEN
          WRITE (numout,'(A39)',advance='no') &
@@ -688 +730 @@
      END IF
 
+#endif
+
    END SUBROUTINE partition_mca_rk
 
@@ -693 +737 @@
    SUBROUTINE partition_rk_core( mask, nx, ny, maxfax, fx, nfx, fy, nfy,   &
                                  ierr )
+#if defined key_mpp_mpi
        USE lib_mpp, ONLY: mppsize
+#endif
        IMPLICIT NONE
        !!------------------------------------------------------------------
@@ -1171 +1217 @@
 
                  IF ( depth(i,j) == 1 ) THEN
-                     newbound = MAX(i - jpreci, pielb(iproc))
+                     newbound = MAX(i - jpreci - nextra, pielb(iproc))
 #if defined TRIM_DEBUG
                  IF ( lwp ) THEN
@@ -1213 +1259 @@
                      ! We've found a wet point in this column so this is as far 
                      ! as we can trim.
-                     newbound = MIN(i + jpreci, pieub(iproc))
+                     newbound = MIN(i + jpreci + nextra, pieub(iproc))
 #if defined TRIM_DEBUG
                      IF ( lwp ) THEN
@@ -1255 +1301 @@
              DO i=MAX(1,pielb(iproc)-jpreci),MIN(jpiglo,pieub(iproc)+jpreci)
                 IF ( depth(i,j) == 1) THEN
-                   newbound = MAX(j - jpreci, pjelb(iproc))
+                   newbound = MAX(j - jpreci - nextra, pjelb(iproc))
 #if defined TRIM_DEBUG
                    IF ( lwp ) THEN
@@ -1298 +1344 @@
             DO i=MAX(1,pielb(iproc)-jpreci),MIN(jpiglo,pieub(iproc)+jpreci)
                IF ( depth(i,j) == 1 ) THEN
-                  newbound = MIN(j + jpreci, pjeub(iproc))
+                  newbound = MIN(j + jpreci + nextra, pjeub(iproc))
 #if defined TRIM_DEBUG
                   IF ( lwp ) then
@@ -1371 +1417 @@
         END IF
 
+#if defined key_mpp_mpi
         IF ( nn_pttrim ) THEN
+           nextra = 2
            CALL part_trim ( imask, trimmed, ierr )
         ELSE
+           ! Need non-zero nextra because otherwise hit trouble with fields
+           ! not being read from disk over land regions
+           nextra = 2
+           !nextra = 0 ! Don't need to back-off on message trimming
+                      ! if we're not trimming the domains
            trimmed(1:4,1:nprocp) = .FALSE.
         ENDIF
+#else
+        trimmed(1:4,1:nprocp) = .FALSE.
+#endif
 
         ! Lower boundary (long.) of sub-domain, GLOBAL coords
@@ -1475 +1531 @@
            njmpp = njmpp - jprecj
         END IF
+        ! ARPDBG - should we allow for trimming of northern edge of
+        ! sub-domains here?
         jubext = pjubext(narea)
         IF(jubext)THEN
@@ -1481 +1539 @@
         END IF
 
-      jelb   = pjelb (narea) ! Lower bound of internal domain
-      jeub   = pjeub (narea) ! Upper bound of internal domain
-      jesub  = pjesub(narea) ! Extent of internal domain
-
-      jpj  = jesub + 2*jprecj ! jpj is the same for all domains - this is
-                              ! what original decomposition did
-      nlcj = jpj
+        jelb   = pjelb (narea) ! Lower bound of internal domain
+        jeub   = pjeub (narea) ! Upper bound of internal domain
+        jesub  = pjesub(narea) ! Extent of internal domain
+
+        jpj  = jesub + 2*jprecj ! jpj is the same for all domains - this is
+                                ! what original decomposition did
+        nlcj = jpj
 
 ! Unlike the East-West boundaries, the global domain does not include
@@ -1540 +1598 @@
    END SUBROUTINE finish_partition
 
-!!$ ARPDBG - we don't want to change the North-fold code for the minute
-!!$     SUBROUTINE mpp_ini_north
-!!$    !!----------------------------------------------------------------------
-!!$    !!               ***  routine mpp_ini_north  ***
-!!$    !!
-!!$    !! ** Purpose :   Initialize special communicator for north folding 
-!!$    !!      condition together with global variables needed in the mpp folding
-!!$    !!
-!!$    !! ** Method  : - Look for northern processors
-!!$    !!              - Put their number in nrank_north
-!!$    !!              - Create groups for the world processors and the north processors
-!!$    !!              - Create a communicator for northern processors
-!!$    !!
-!!$    !! ** output
-!!$    !!      njmppmax = njmpp for northern procs
-!!$    !!      ndim_rank_north = number of processors in the northern line
-!!$    !!      nrank_north (ndim_rank_north) = number  of the northern procs.
-!!$    !!      ngrp_world = group ID for the world processors
-!!$    !!      ngrp_north = group ID for the northern processors
-!!$    !!      ncomm_north = communicator for the northern procs.
-!!$    !!      north_root = number (in the world) of proc 0 in the northern comm.
-!!$    !!      nwidthmax = width of widest northern domain
-!!$    !!
-!!$    !! History :
-!!$    !!        !  03-09 (J.M. Molines, MPI only )
-!!$    !!        !  08-09 (A.R. Porter - for new decomposition)
-!!$    !!----------------------------------------------------------------------
-!!$    USE exchmod, ONLY: nrank_north, north_root, ndim_rank_north, &
-!!$                       ncomm_north, ngrp_world, ngrp_north
-!!$    IMPLICIT none
-!!$#ifdef key_mpp_shmem
-!!$    CALL ctl_stop( ' mpp_ini_north not available in SHMEM' )
-!!$# elif key_mpp_mpi
-!!$    INTEGER :: ierr
-!!$    INTEGER :: jproc
-!!$    INTEGER :: ii,ji
-!!$    !!----------------------------------------------------------------------
-!!$
-!!$    ! Look for how many procs on the northern boundary
-!!$    !
-!!$    ndim_rank_north = 0
-!!$    nwidthmax = 0
-!!$
-!!$    DO jproc=1,jpnij
-!!$       IF ( pjubext(jproc) ) THEN
-!!$          ndim_rank_north = ndim_rank_north + 1
-!!$
-!!$          ! and for the width of the widest northern domain...
-!!$          IF(piesub(jproc) > nwidthmax)THEN
-!!$             nwidthmax = piesub(jproc)
-!!$          END IF
-!!$       END IF
-!!$    END DO
-!!$    nwidthmax = nwidthmax + 2*jpreci ! Allow for halos
-!!$
-!!$    ! Allocate the right size to nrank_north
-!!$    !
-!!$    ALLOCATE(nrank_north(ndim_rank_north))
-!!$
-!!$    ! Fill the nrank_north array with proc. number of northern procs.
-!!$    ! Note : the rank start at 0 in MPI
-!!$    !
-!!$    ii=0
-!!$    DO ji = 1, jpnij
-!!$       IF ( pjubext(ji) ) THEN
-!!$          ii=ii+1
-!!$          nrank_north(ii)=ji-1
-!!$       END IF
-!!$    END DO
-!!$    ! create the world group
-!!$    !
-!!$    CALL MPI_COMM_GROUP(mpi_comm_opa,ngrp_world,ierr)
-!!$    !
-!!$    ! Create the North group from the world group
-!!$    CALL MPI_GROUP_INCL(ngrp_world,ndim_rank_north,nrank_north,ngrp_north,ierr)
-!!$
-!!$    ! Create the North communicator , ie the pool of procs in the north group
-!!$    !
-!!$    CALL MPI_COMM_CREATE(mpi_comm_opa,ngrp_north,ncomm_north,ierr)
-!!$
-!!$
-!!$    ! find proc number in the world of proc 0 in the north
-!!$    CALL MPI_GROUP_TRANSLATE_RANKS(ngrp_north,1,0,ngrp_world,north_root,ierr)
-!!$#endif
-!!$
-!!$  END SUBROUTINE mpp_ini_north
-
-     SUBROUTINE eval_partition( nx, ny, mask, score )
-
-        ! Compute the cost function for the current partition
-        !
-        ! Assume that the time taken for a run is proportional
-        ! to the maximum over processors of:
-        !     w_processing * cost_processing
-        !   + w_communications * cost_communications
-        ! Assume further that cost_processing goes as
-        !   (number of wet points) + f_proc * (number of dry points)
-        ! (with f_proc << 1)
-        ! and that cost_communications goes as
-        !   (cost of intra-node communications) +
-        !   f_comm * (cost of inter-node communications)
-        ! (with f_comm << 1)
-        !
-        ! However, because of the possiblity of network contention,
-        ! other factors may also matter, especially:
-        !   total over sub-domains of halo points with off-node neighbours
-        !   max over nodes of total off-node halo points and message counts
-        !
-        ! With this in mind, we construct the ansatz
-        !  maximum over processors of {
-        !     w_1 * (number of wet points)
-        !   + w_2 * (number of dry points)
-        !   + w_3 * (halo points with off-node neighbours)
-        !   + w_4 * (halo points with on-node neighbours)
-        !   + ...
-        ! }
-        USE lib_mpp,     ONLY: mppsize
-        USE mapcomm_mod, ONLY: iprocmap, land
-        IMPLICIT NONE
-        !     Arguments
+
+   SUBROUTINE mpp_ini_north
+      !!----------------------------------------------------------------------
+      !!               ***  routine mpp_ini_north  ***
+      !!
+      !! ** Purpose :   Initialize special communicator for north folding 
+      !!      condition together with global variables needed in the mpp folding
+      !!
+      !! ** Method  : - Look for northern processors
+      !!              - Put their number in nrank_north
+      !!              - Create groups for the world processors and the north 
+      !!                processors
+      !!              - Create a communicator for northern processors
+      !!
+      !! ** output
+      !!      njmppmax = njmpp for northern procs
+      !!      ndim_rank_north = number of processors in the northern line
+      !!      nrank_north (ndim_rank_north) = number  of the northern procs.
+      !!      ngrp_world = group ID for the world processors
+      !!      ngrp_north = group ID for the northern processors
+      !!      ncomm_north = communicator for the northern procs.
+      !!      north_root = number (in the world) of proc 0 in the northern comm.
+      !!      nwidthmax = width of widest northern domain
+      !!
+      !! History :
+      !!        !  03-09 (J.M. Molines, MPI only )
+      !!        !  08-09 (A.R. Porter - for new decomposition)
+      !!----------------------------------------------------------------------
+      USE par_oce, ONLY: jperio, jpni
+      USE exchmod, ONLY: nrank_north, north_root, ndim_rank_north, &
+                         ncomm_north, ngrp_world, ngrp_north,      &
+                         do_nfold, num_nfold_rows, nfold_npts
+      USE dom_oce, ONLY: narea
+      IMPLICIT none
+#ifdef key_mpp_shmem
+      CALL ctl_stop('STOP', ' mpp_ini_north not available in SHMEM' )
+# elif key_mpp_mpi
+      INTEGER :: ierr
+      INTEGER :: jproc
+      INTEGER :: ii,ji
+      !!----------------------------------------------------------------------
+
+      ! Look for how many procs on the northern boundary
+      !
+      ndim_rank_north = 0
+      nwidthmax       = 0
+      do_nfold        = .FALSE.
+
+      IF (.NOT. (jperio >= 3 .AND. jperio <= 6 .AND. jpni > 1) ) THEN
+         ! No northern boundary to worry about
+         RETURN
+      END IF
+
+      DO jproc=1,mppsize,1
+         IF ( pjubext(jproc) ) THEN
+
+            ! If trimming of dry land from sub-domains is enabled
+            ! then check that this PE does actually have data to
+            ! contribute to the N-fold. If trimming is not enabled
+            ! then this condition will always be true for northern
+            ! PEs.
+            IF( pjeub(jproc) > (jpjglo - num_nfold_rows) )THEN
+
+               ndim_rank_north = ndim_rank_north + 1
+
+               ! and for the width of the widest northern domain...
+               nwidthmax = MAX(nwidthmax, piesub(jproc))
+            ENDIF
+
+         END IF
+      END DO
+      nwidthmax = nwidthmax + 2*jpreci ! Allow for halos
+
+      ! Allocate the right size to nrank_north
+      !
+      ALLOCATE(nrank_north(ndim_rank_north), nfold_npts(ndim_rank_north), &
+               Stat=ierr)
+      IF( ierr /= 0 )THEN
+         CALL ctl_stop('STOP','mpp_ini_north: failed to allocate arrays')
+      END IF
+
+#if defined PARTIT_DEBUG
+      IF(lwp)THEN
+         WRITE(*,*) 'mpp_ini_north: no. of northern PEs = ',ndim_rank_north
+         WRITE(*,*) 'mpp_ini_north: nwidthmax = ',nwidthmax
+      END IF
+#endif
+      ! Fill the nrank_north array with proc. number of northern procs.
+      ! Note : ranks start at 0 in MPI
+      !
+      ii=0
+      DO ji = 1, mppsize, 1
+         IF (  pjubext(ji)       .AND.          &
+              (pjeub(ji) > (jpjglo - num_nfold_rows)) ) THEN
+            ii=ii+1
+            nrank_north(ii)=ji-1
+
+            ! Flag that this PE does do North-fold (with trimming, checking
+            ! npolj is no longer sufficient)
+            IF(ji == narea) do_nfold = .TRUE.
+
+#if defined NO_NFOLD_GATHER
+            ! How many data points will this PE have to send for N-fold?
+
+            ! No. of valid rows for n-fold = num_nfold_rows - <no. trimmed rows>
+            !                              = num_nfold_rows - jpjglo + pjeub(ji)
+
+            ! ARPDBG - could trim land-only rows/cols from this...
+            nfold_npts(ii) = MAX(num_nfold_rows - jpjglo + pjeub(ji), 0) * &
+                             ( nleit(ji) - nldit(ji) + 1 )
+#endif 
+         END IF
+      END DO
+      ! create the world group
+      !
+      CALL MPI_COMM_GROUP(mpi_comm_opa,ngrp_world,ierr)
+      !
+      ! Create the North group from the world group
+      CALL MPI_GROUP_INCL(ngrp_world,ndim_rank_north,nrank_north, &
+                          ngrp_north,ierr)
+
+      ! Create the North communicator , ie the pool of procs in the north group
+      !
+      CALL MPI_COMM_CREATE(mpi_comm_opa,ngrp_north,ncomm_north,ierr)
+
+
+      ! find proc number in the world of proc 0 in the north
+      CALL MPI_GROUP_TRANSLATE_RANKS(ngrp_north,1,0,ngrp_world,north_root,ierr)
+
+#endif
+
+   END SUBROUTINE mpp_ini_north
+
+
+   SUBROUTINE eval_partition( nx, ny, mask, score )
+
+      ! Compute the cost function for the current partition
+      !
+      ! Assume that the time taken for a run is proportional
+      ! to the maximum over processors of:
+      !     w_processing * cost_processing
+      !   + w_communications * cost_communications
+      ! Assume further that cost_processing goes as
+      !   (number of wet points) + f_proc * (number of dry points)
+      ! (with f_proc << 1)
+      ! and that cost_communications goes as
+      !   (cost of intra-node communications) +
+      !   f_comm * (cost of inter-node communications)
+      ! (with f_comm << 1)
+      !
+      ! However, because of the possiblity of network contention,
+      ! other factors may also matter, especially:
+      !   total over sub-domains of halo points with off-node neighbours
+      !   max over nodes of total off-node halo points and message counts
+      !
+      ! With this in mind, we construct the ansatz
+      !  maximum over processors of {
+      !     w_1 * (number of wet points)
+      !   + w_2 * (number of dry points)
+      !   + w_3 * (halo points with off-node neighbours)
+      !   + w_4 * (halo points with on-node neighbours)
+      !   + ...
+      ! }
+#if defined key_mpp_mpi
+      USE lib_mpp,     ONLY: mppsize
+#endif
+      USE mapcomm_mod, ONLY: iprocmap, land
+      IMPLICIT NONE
+      !     Arguments
         INTEGER, INTENT(in) :: nx, ny
         INTEGER, INTENT(in) :: mask(nx,ny)
@@ -1704 +1811 @@
         ! next nn_cpnode ranks are assigned to node 1, etc
         INTEGER, ALLOCATABLE :: node(:)
+
+#if defined key_mpp_mpi
 
         ALLOCATE(node(nprocp))
@@ -1940 +2049 @@
       DEALLOCATE(node)
 
+#endif
+
      END SUBROUTINE eval_partition
 
@@ -2395 +2506 @@
 
 
-    SUBROUTINE smooth_bathy(inbathy)
+    SUBROUTINE smooth_global_bathy(inbathy, imask)
        USE dom_oce
-       USE domzgr
+       USE domzgr, ONLY: rn_sbot_min, rn_sbot_max, rn_theta, rn_thetb, &
+                         rn_rmax, ln_s_sigma, rn_bb, rn_hc, fssig1, &
+                         namzgr_sco
        USE in_out_manager, ONLY: numnam
        IMPLICIT none
        !!----------------------------------------------------------------------
-       !!                      Routine smooth_bathy
+       !!                      Routine smooth_global_bathy
        !!   Replicates the smoothing done on the decomposed domain in zgr_sco()
        !!   in domzgr.F90. However, here the domain is NOT decomposed and
@@ -2408 +2521 @@
        !!   is done using a mask that is the same as that which is eventually
        !!   computed after zgr_sco() has been called. (The smoothing process
-       !!   below can change whether grid points are wet or dry.)
+       !!   below can (erroneously) change whether grid points are wet or dry.)
        !!----------------------------------------------------------------------
        REAL(wp), INTENT(inout), DIMENSION(:,:) :: inbathy ! The bathymetry to 
                                                           ! be smoothed
+       INTEGER, INTENT(inout), DIMENSION(:,:)  :: imask   ! Mask holding index of
+                                                          ! bottom level
        ! Locals
-       INTEGER  :: ji, jj, jl, ierr
+       INTEGER  :: ji, jj, jk, jl, ierr
        INTEGER  :: iip1, ijp1, iim1, ijm1   ! temporary integers
        INTEGER  :: x_size, y_size
-       REAL(wp) :: zrmax, zri, zrj
+       REAL(wp) :: zrmax, zri, zrj, zcoeft
        REAL(wp), PARAMETER :: TOL_ZERO = 1.0E-20_wp ! Any value less than 
                                                     ! this assumed zero
-       REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zenv, ztmp, zmsk, zbot
+       REAL(wp), ALLOCATABLE, DIMENSION(:,:)   :: zenv, ztmp, zmsk, zbot, &
+                                                  zscosrf, zhbatt
+       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zgsigt3, zgdept
        !
-       NAMELIST/namzgr_sco/ rn_sbot_max, rn_sbot_min, rn_theta, rn_thetb, &
-                            rn_rmax, ln_s_sigma, rn_bb, rn_hc
        !!----------------------------------------------------------------------
 
@@ -2431 +2546 @@
 
        ALLOCATE(zenv(x_size,y_size), ztmp(x_size,y_size), zmsk(x_size,y_size), &
-                zbot(x_size,y_size), Stat=ierr)
+                zbot(x_size,y_size), zgdept(x_size,y_size,jpkdta), zhbatt(x_size, y_size), &
+                zscosrf(x_size,y_size), zgsigt3(x_size,y_size,jpkdta), Stat=ierr)
        IF( ierr /= 0 ) THEN
-          CALL ctl_stop('smooth_bathy: ERROR - failed to allocate workspace arrays')
+          CALL ctl_stop('smooth_global_bathy: ERROR - failed to allocate workspace arrays')
           RETURN
        ENDIF
@@ -2441 +2557 @@
        READ  ( numnam, namzgr_sco )
 
+       zscosrf(:,:) = 0._wp            ! ocean surface depth (here zero: no under ice-shelf sea)
        zbot(:,:) = inbathy(:,:)        ! ocean bottom depth
        !                               ! set maximum ocean depth
@@ -2477 +2594 @@
 
          !
-         IF(lwp)WRITE(numout,"('smooth_bathy : iter=',I5,' rmax=',F8.4,' nb of pt= ',I8)") &
+         IF(lwp)WRITE(numout,"('smooth_global_bathy : iter=',I5,' rmax=',F8.4,' nb of pt= ',I8)") &
                                                          jl, zrmax, INT( SUM(zmsk(:,:) ) )
          !
@@ -2491 +2608 @@
                ijm1 = MAX( jj-1,  1  )      ! first raw  (jj=nlcj)
                IF( zmsk(ji,jj) == 1._wp ) THEN
-                  ztmp(ji,jj) =   (                                                                                   &
-             &      zenv(iim1,ijp1)*zmsk(iim1,ijp1) + zenv(ji,ijp1)*zmsk(ji,ijp1) + zenv(iip1,ijp1)*zmsk(iip1,ijp1)   &
-             &    + zenv(iim1,jj  )*zmsk(iim1,jj  ) + zenv(ji,jj  )*    2._wp     + zenv(iip1,jj  )*zmsk(iip1,jj  )   &
-             &    + zenv(iim1,ijm1)*zmsk(iim1,ijm1) + zenv(ji,ijm1)*zmsk(ji,ijm1) + zenv(iip1,ijm1)*zmsk(iip1,ijm1)   &
-             &                    ) / (                                                                               &
-             &                      zmsk(iim1,ijp1) +               zmsk(ji,ijp1) +                 zmsk(iip1,ijp1)   &
-             &    +                 zmsk(iim1,jj  ) +                   2._wp     +                 zmsk(iip1,jj  )   &
-             &    +                 zmsk(iim1,ijm1) +               zmsk(ji,ijm1) +                 zmsk(iip1,ijm1)   &
-             &                        )
+                ztmp(ji,jj) =   (                                                                                  &
+             &    zenv(iim1,ijp1)*zmsk(iim1,ijp1) + zenv(ji,ijp1)*zmsk(ji,ijp1) + zenv(iip1,ijp1)*zmsk(iip1,ijp1)  &
+             &  + zenv(iim1,jj  )*zmsk(iim1,jj  ) + zenv(ji,jj  )*    2._wp     + zenv(iip1,jj  )*zmsk(iip1,jj  )  &
+             &  + zenv(iim1,ijm1)*zmsk(iim1,ijm1) + zenv(ji,ijm1)*zmsk(ji,ijm1) + zenv(iip1,ijm1)*zmsk(iip1,ijm1)  &
+             &                  ) / (                                                                              &
+             &                    zmsk(iim1,ijp1) +               zmsk(ji,ijp1) +                 zmsk(iip1,ijp1)  &
+             &  +                 zmsk(iim1,jj  ) +                   2._wp     +                 zmsk(iip1,jj  )  &
+             &  +                 zmsk(iim1,ijm1) +               zmsk(ji,ijm1) +                 zmsk(iip1,ijm1)  &
+             &                      )
                ENDIF
             END DO
@@ -2514 +2631 @@
       !                                                   ! ================ !
       !
+      !                                        ! envelop bathymetry saved in zhbatt
+      zhbatt(:,:) = zenv(:,:) 
+      ! gphit calculated in nemo_init->dom_init->dom_hgr and dom_hgr requires that 
+      ! partitioning already done. Could repeat its calculation here but since AMM doesn't
+      ! require it we leave it out for the moment ARPDBG
+      CALL ctl_warn( ' ARPDBG - NOT checking whether s-coordinates are tapered in vicinity of the Equator' )
+!!$      IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN
+!!$         CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' )
+!!$         DO jj = 1, jpj
+!!$            DO ji = 1, jpi
+!!$               ztaper = EXP( -(gphit(ji,jj)/8._wp)**2 )
+!!$               hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1._wp - ztaper )
+!!$            END DO
+!!$         END DO
+!!$      ENDIF
+
       ! Subtract off rn_sbot_min so can check for land using zenv = LAND (0)
       inbathy(:,:) = zenv(:,:) - rn_sbot_min
 
-!!$      IF(lwp)THEN
-!!$         OPEN(UNIT=1098, FILE='bathy_bottom.dat', STATUS='REPLACE', &
-!!$              ACTION='WRITE', IOSTAT=jj)
-!!$         IF(jj == 0)THEN
-!!$            DO jj = 1, y_size
-!!$               DO ji = 1, x_size
-!!$                  WRITE (1098,"(I4,1x,I4,3(E14.4,1x))") ji, jj, &
-!!$                        inbathy(ji,jj), zbot(ji,jj), &
-!!$                       (inbathy(ji,jj)-zbot(ji,jj))
-!!$               END DO
-!!$               WRITE (1098,*)
-!!$            END DO
-!!$            CLOSE(1098)
-!!$         END IF
-!!$      END IF
-
-    END SUBROUTINE smooth_bathy
+
+      !                                            ! =======================
+      !                                            !   s-ccordinate fields     (gdep., e3.)
+      !                                            ! =======================
+      !
+      ! non-dimensional "sigma" for model level depth at w- and t-levels
+
+      IF( ln_s_sigma ) THEN        ! Song and Haidvogel style stretched sigma for depths
+         !                         ! below rn_hc, with uniform sigma in shallower waters
+         DO ji = 1, x_size
+            DO jj = 1, y_size
+
+               IF( zhbatt(ji,jj) > rn_hc ) THEN    !deep water, stretched sigma
+                  DO jk = 1, jpk
+                     zgsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp)       , rn_bb )
+                  END DO
+               ELSE ! shallow water, uniform sigma
+                  DO jk = 1, jpk
+                     zgsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp)
+                  END DO
+               ENDIF
+               !
+               DO jk = 1, jpk
+                  zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
+                  zgdept (ji,jj,jk) = zscosrf(ji,jj) + (zhbatt(ji,jj)-rn_hc)*zgsigt3(ji,jj,jk)+rn_hc*zcoeft
+               END DO
+               !
+            END DO   ! for all jj's
+         END DO    ! for all ji's
+      ELSE
+         CALL ctl_stop('STOP', &
+                       'partition_mod::smooth_global_bathy() only supports ln_s_sigma = .TRUE. currently!')
+      END IF
+
+      ! HYBRID scheme
+      DO jj = 1, y_size
+         DO ji = 1, x_size
+            DO jk = 1, jpkm1
+               IF( zbot(ji,jj) >= zgdept(ji,jj,jk) )  imask(ji,jj) = MAX( 2, jk )
+               IF( zbot(ji,jj) == 0._wp           )   imask(ji,jj) = 0
+            END DO
+         END DO
+      END DO
+
+      ! Dump to file for debugging ARPDBG
+      IF(lwp)THEN
+         OPEN(UNIT=1098, FILE='smoothed_bathy.dat', STATUS='REPLACE', &
+              ACTION='WRITE', IOSTAT=jj)
+         IF(jj == 0)THEN
+            DO jj = 1, y_size
+               DO ji = 1, x_size
+                  WRITE (1098,"(I4,1x,I4,3(E14.4,1x),I4)") ji, jj, &
+                        inbathy(ji,jj),             zbot(ji,jj),   &
+                        inbathy(ji,jj)-zbot(ji,jj), imask(ji,jj)
+               END DO
+               WRITE (1098,*)
+            END DO
+            CLOSE(1098)
+         END IF
+      END IF
+
+    END SUBROUTINE smooth_global_bathy
+
+
+    SUBROUTINE global_bot_level(imask)
+      USE par_oce, ONLY: jperio
+      IMPLICIT none
+      !!----------------------------------------------------------------------
+      !! Compute the deepest level for any of the u,v,w or T grids. (Code
+      !! taken from zgr_bot_level() and intermediate arrays for U and V
+      !! removed.)
+      !!----------------------------------------------------------------------
+      INTEGER, DIMENSION(:,:), INTENT(inout) :: imask
+      ! Locals
+      INTEGER :: ji, jj
+      INTEGER :: x_size, y_size
+
+       ! Do this because we've not decomposed the domain yet and therefore
+       ! jpi,jpj,nlc{i,j} etc. are not set.
+       x_size = SIZE(imask, 1)
+       y_size = SIZE(imask, 2)
+
+      imask(:,:) = MAX( imask(:,:) , 1 )  ! bottom k-index of T-level (=1 over land)
+
+      !
+      ! Compute and store the deepest bottom k-index of any grid-type at 
+      ! each grid point.
+      ! For use in removing data below ocean floor from halo exchanges.
+      DO jj = 1, y_size-1
+         DO ji = 1, x_size-1
+            imask(ji,jj) = MAX(imask(ji,jj)+1,                           & ! W (= T-level + 1)
+                               MIN(  imask(ji+1,jj  ) , imask(ji,jj)  ), & ! U
+                               MIN(  imask(ji  ,jj+1) , imask(ji,jj)  ) )  ! V
+         END DO
+         imask(x_size,jj) = imask(x_size-1,jj)
+      END DO
+
+      ! Check on jperio because we've not set cyclic_bc in mapcomms yet
+      IF(jperio == 1 .OR. jperio == 4 .OR. jperio == 6)THEN
+         ! Impose global cyclic boundary conditions on the array holding the
+         ! deepest level
+         imask(1,:)      = imask(x_size - 1, :)
+         imask(x_size,:) = imask(2,:)
+      END IF
+
+      ! Dump to file for debugging ARPDBG
+      IF(lwp)THEN
+         OPEN(UNIT=1098, FILE='bathy_bottom.dat', STATUS='REPLACE', &
+              ACTION='WRITE', IOSTAT=jj)
+         IF(jj == 0)THEN
+            DO jj = 1, y_size
+               DO ji = 1, x_size
+                  WRITE (1098,"(I4,1x,I4,1x,I4)") ji, jj, imask(ji,jj)
+               END DO
+               WRITE (1098,*)
+            END DO
+            CLOSE(1098)
+         END IF
+      END IF
+
+    END SUBROUTINE global_bot_level
 
 END MODULE partition_mod

branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn_c3d.h90

@@ -281 +281 @@
          zahmeq = 5.0 * aht0
          zahmm  = min( 160000.0, ahm0)
-         zemax = MAXVAL ( e1t(:,:) * e2t(:,:), tmask(:,:,1) .GE. 0.5 )
-         zemin = MINVAL ( e1t(:,:) * e2t(:,:), tmask(:,:,1) .GE. 0.5 )
-         zeref = MAXVAL ( e1t(:,:) * e2t(:,:),   &
-             &   tmask(:,:,1) .GE. 0.5 .AND. ABS(gphit(:,:)) .GT. 50. )
+         zemax = MAXVAL ( e1t(:,:) * e2t(:,:), MASK=(tmask(:,:,1) .GE. 0.5) )
+         zemin = MINVAL ( e1t(:,:) * e2t(:,:), MASK=(tmask(:,:,1) .GE. 0.5) )
+         zeref = MAXVAL ( e1t(:,:) * e2t(:,:),                 &
+                          MASK=( (tmask(:,:,1) .GE. 0.5) .AND. &
+                                 (ABS(gphit(:,:)) .GT. 50.) ) )
  
          DO jj = 1, jpj

branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90

@@ -706 +706 @@
          IF( LEN( TRIM(sdf_n(jf)%wname) ) > 0 )   sdf(jf)%wgtname = TRIM( cdir )//TRIM( sdf_n(jf)%wname )
          sdf(jf)%vcomp   = sdf_n(jf)%vcomp
+         sdf(jf)%rotn    = .FALSE.
+         ! Initialise arrays allocated in calling routine
+         sdf(jf)%fnow(:,:,:)                              = 0.0_wp
+         IF(ALLOCATED(sdf(jf)%fdta))sdf(jf)%fdta(:,:,:,:) = 0.0_wp
       END DO
 

branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90

@@ -179 +179 @@
          END DO
          !                                         ! fill sf with slf_i and control print
-         CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_core', 'flux formulation for ocean surface boundary condition', 'namsbc_core' )
+         CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_core', &
+                        'flux formulation for ocean surface boundary condition', &
+                        'namsbc_core' )
          !
       ENDIF
@@ -258 +260 @@
          END DO
       END DO
+
       CALL lbc_lnk( zwnd_i(:,:) , 'T', -1. )
       CALL lbc_lnk( zwnd_j(:,:) , 'T', -1. )
@@ -346 +349 @@
          END DO
       END DO
+
       CALL lbc_lnk( utau(:,:), 'U', -1. )
       CALL lbc_lnk( vtau(:,:), 'V', -1. )

branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/SOL/solpcg.F90

@@ -89 +89 @@
       USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
       USE wrk_nemo, ONLY:   zgcr => wrk_2d_1
+!      USE arpdebugging, ONLY: dump_array
       !!
       INTEGER, INTENT(inout) ::   kindic   ! solver indicator, < 0 if the conver-
@@ -97 +98 @@
       REAL(wp) ::   zgcad        ! temporary scalars
       REAL(wp), DIMENSION(2) ::   zsum
+      INTEGER, SAVE :: istep = 0 ! ARPDBG
       !!----------------------------------------------------------------------
       
@@ -107 +109 @@
       zgcr = 0._wp
       gcr  = 0._wp
+!      CALL dump_array(istep, 'gcx_pre_lbc', gcx, withHalos=.TRUE.)
 
       CALL lbc_lnk( gcx, c_solver_pt, 1. )   ! lateral boundary condition
+
+      istep = istep + 1
+!      CALL dump_array(istep, 'gcx', gcx, withHalos=.TRUE.)
+!      CALL dump_array(istep, 'gcp', gcp(:,:,1), withHalos=.TRUE.)
+!      CALL dump_array(istep, 'gcb', gcb, withHalos=.TRUE.)
+!      CALL dump_array(istep, 'ua', ua(:,:,1), withHalos=.TRUE.)
 
       ! gcr   = gcb-a.gcx

branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90

@@ -98 +98 @@
       !! ** Action :   Update pta arrays with the before rotated diffusion
       !!----------------------------------------------------------------------
+!      USE arpdebugging, ONLY: dump_array
       USE timing,   ONLY: timing_start, timing_stop
       USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
@@ -158 +159 @@
       ENDIF
       !
+!      CALL dump_array(kt, 'ptb', ptb(:,:,1,1), withHalos=.TRUE.)
       !                                                          ! ===========
 !DIR$ SHORTLOOP

branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90

@@ -235 +235 @@
       IF( zdf_ddm_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_ddm_init : unable to allocate arrays' )
       !
+      ! Initialise new array
+      avs(:,:,:) = 0.0_wp
+
    END SUBROUTINE zdf_ddm_init
 

branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/arpdebugging.f90

@@ -11 +11 @@
   CONTAINS
 
-    SUBROUTINE dump_rarray(count, name, field1, field2, withHalos)
+    SUBROUTINE dump_rarray(count, name, field1, field2, withHalos, &
+                           toGlobal)
       IMPLICIT none
       INTEGER,           INTENT(in) :: count  ! What timestep we're on
@@ -18 +19 @@
       REAL,              INTENT(in), DIMENSION(:,:), OPTIONAL :: field2
       LOGICAL,           INTENT(in),                 OPTIONAL :: withHalos
+      LOGICAL,           INTENT(in),                 OPTIONAL :: toGlobal
       ! Locals
       INTEGER           :: ji, jj
       CHARACTER (len=4) :: crank,ccount
-      LOGICAL           :: lwithHalos
+      LOGICAL           :: lwithHalos, ltoGlobal
       INTEGER           :: ibound, jbound
       INTEGER, DIMENSION(2) :: shape1, shape2
@@ -37 +39 @@
       lwithHalos = .false.
       IF(present(withHalos))lwithHalos = withHalos
+      ! By default we convert to global coordinates rather than those local
+      ! to this process
+      ltoGlobal = .true.
+      IF(present(toGlobal))ltoGlobal = toGlobal
 
       WRITE(crank,FMT="(I4)") narea-1
@@ -57 +63 @@
             END DO
          ELSE
-            DO jj=nldj,nlej,1
-               DO ji=nldi,nlei,1
-                  WRITE(997,FMT=fmt_var) mig(ji),mjg(jj), &
-                                               field1(ji,jj)
-               END DO
-            END DO
+
+            IF(ltoGlobal)THEN
+               DO jj=nldj,nlej,1
+                  DO ji=nldi,nlei,1
+                     WRITE(997,FMT=fmt_var) mig(ji),mjg(jj), &
+                                            field1(ji,jj)
+                  END DO
+               END DO
+            ELSE
+               DO jj=nldj,nlej,1
+                  DO ji=nldi,nlei,1
+                     WRITE(997,FMT=fmt_var) ji,jj, field1(ji,jj)
+                  END DO
+               END DO
+            END IF
          END IF
 
@@ -71 +86 @@
 
          IF(PRESENT(field2))THEN
-            DO ji=1,ibound,1
-               DO jj=1,jbound,1
-                  WRITE(997,FMT=fmt_var) mig(ji), mjg(jj), &
-                                                field1(ji,jj),field2(ji,jj)
-               END DO
-               WRITE(997,*)
-            END DO
-         ELSE
-            DO ji=1,ibound,1
-               DO jj=1,jbound,1
-                  WRITE(997,FMT=fmt_var) mig(ji), mjg(jj), field1(ji,jj)
-               END DO
-               WRITE(997,*)
-            END DO
+            IF(ltoGlobal)THEN
+               DO ji=1,ibound,1
+                  DO jj=1,jbound,1
+                     WRITE(997,FMT=fmt_var) mig(ji), mjg(jj), &
+                                            field1(ji,jj),field2(ji,jj)
+                  END DO
+                  WRITE(997,*)
+               END DO
+            ELSE
+               DO ji=1,ibound,1
+                  DO jj=1,jbound,1
+                     WRITE(997,FMT=fmt_var) ji, jj, field1(ji,jj), field2(ji,jj)
+                  END DO
+                  WRITE(997,*)
+               END DO
+            END IF
+         ELSE
+            IF(ltoGlobal)THEN
+               DO ji=1,ibound,1
+                  DO jj=1,jbound,1
+                     WRITE(997,FMT=fmt_var) mig(ji), mjg(jj), field1(ji,jj)
+                  END DO
+                  WRITE(997,*)
+               END DO
+            ELSE
+               DO ji=1,ibound,1
+                  DO jj=1,jbound,1
+                     WRITE(997,FMT=fmt_var) ji, jj, field1(ji,jj)
+                  END DO
+                  WRITE(997,*)
+               END DO
+            END IF
          END IF
 
@@ -94 +127 @@
 
 
-    SUBROUTINE dump_iarray(count, name, field1, field2, withHalos)
+    SUBROUTINE dump_iarray(count, name, field1, field2, withHalos, &
+                           toGlobal)
       IMPLICIT none
       INTEGER,           INTENT(in) :: count  ! What timestep we're on
@@ -101 +135 @@
       INTEGER,           INTENT(in), DIMENSION(:,:), OPTIONAL :: field2
       LOGICAL,           INTENT(in),                 OPTIONAL :: withHalos
+      LOGICAL,           INTENT(in),                 OPTIONAL :: toGlobal
       ! Locals
       INTEGER           :: ji, jj
       CHARACTER (len=4) :: crank,ccount
-      LOGICAL           :: lwithHalos
+      LOGICAL           :: lwithHalos, ltoGlobal
       INTEGER           :: ibound, jbound
       INTEGER, DIMENSION(2) :: shape1, shape2
@@ -120 +155 @@
       lwithHalos = .false.
       IF(present(withHalos))lwithHalos = withHalos
+      ! By default we convert to global coordinates rather than those local
+      ! to this process
+      ltoGlobal = .true.
+      IF(present(toGlobal))ltoGlobal = toGlobal
 
       WRITE(crank,FMT="(I4)") narea-1
@@ -156 +195 @@
             DO ji=1,ibound,1
                DO jj=1,jbound,1
-                  WRITE(997,FMT=fmt_var) ji, jj, &
+                  WRITE(997,FMT=fmt_var) mig(ji), mjg(jj), &
                                                 field1(ji,jj),field2(ji,jj)
                END DO
@@ -164 +203 @@
             DO ji=1,ibound,1
                DO jj=1,jbound,1
-                  WRITE(997,FMT=fmt_var) ji, jj, field1(ji,jj)
+                  WRITE(997,FMT=fmt_var) mig(ji), mjg(jj), field1(ji,jj)
                END DO
                WRITE(997,*)

branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90

@@ -76 +76 @@
    USE timing, ONLY: timing_init, timing_finalize, timing_disable, timing_enable
 
-!#define ARPDEBUG
+#define ARPDEBUG
 
    IMPLICIT NONE
@@ -235 +235 @@
       lwp = (narea == 1) .OR. ln_ctl                        ! control of all listing output print
 
+      ! Calculate domain z dimensions as needed when partitioning.
+      ! This used to be done in par_oce.F90 when they were parameters rather
+      ! than variables
+      IF( Agrif_Root() ) THEN
+         jpk = jpkdta                                             ! third dim
+         jpkm1 = jpk-1                                            ! inner domain indices
+      ENDIF
+
       CALL timing_init                                      ! Init timing module
       CALL timing_disable                                   ! but disable during startup
@@ -251 +259 @@
          jpnj  = 1
          jpnij = jpni*jpnj
+#endif
+
+#if   defined key_mpp_rkpart
+      ELSE
+         CALL ctl_stop( 'STOP', 'nemo_init : invalid inputs in namelist - cannot specify jpn{i,j}>0 when using recursive k-section paritioning!' )
 #endif
       END IF
@@ -265 +278 @@
          jpij  = jpi*jpj                                          !  jpi x j
 #endif
-         jpk = jpkdta                                             ! third dim
-         jpkm1 = jpk-1                                            ! inner domain indices
       ENDIF
 
@@ -581 +592 @@
 
    SUBROUTINE nemo_recursive_partition( num_pes )
-      USE dom_oce, ONLY: ln_zco, ntopo
-      USE iom,     ONLY: jpiglo, jpjglo, wp, jpdom_unknown, &
-                         iom_open, iom_get, iom_close
+      USE in_out_manager, ONLY: numnam
+      USE dom_oce,        ONLY: ln_zco, ntopo
+      USE dom_oce,        ONLY: gdepw_0, gdept_0, e3w_0, e3t_0, &
+                                mig, mjg, mi0, mi1, mj0, mj1,  mbathy, bathy
+      USE domzgr,         ONLY: zgr_z, zgr_bat, namzgr, zgr_zco, zgr_zps
+      USE closea,         ONLY: dom_clo
+      USE domain,         ONLY: dom_nam
+      USE iom,            ONLY: jpiglo, jpjglo, wp, jpdom_unknown, &
+                                iom_open, iom_get, iom_close
       USE mapcomm_mod, ONLY: ielb, ieub, pielb, pjelb, pieub, pjeub, &
                              iesub, jesub, jeub, ilbext, iubext, jubext, &
                              jlbext, pnactive, piesub, pjesub, jelb, pilbext, &
-                             piubext, pjlbext, pjubext, LAND
-      USE partition_mod, ONLY: partition_rk, partition_mca_rk, imask, smooth_bathy
+                             piubext, pjlbext, pjubext, LAND, msgtrim_z
+      USE partition_mod, ONLY: partition_rk, partition_mca_rk, &
+                               imask, ibotlevel, partition_mask_alloc, &
+                               smooth_global_bathy, global_bot_level
       USE par_oce,       ONLY: do_exchanges
 #if defined key_mpp_mpi
@@ -607 +626 @@
       INTEGER :: ii,jj,iproc                   ! Loop index
       INTEGER :: jparray(2)                    ! Small array for gathering 
+      CHARACTER(LEN=8) :: lstr                 ! Local string for reading env. var.
+      INTEGER          :: lztrim               ! Local int for      "      "    "
       REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdta  ! temporary data workspace
       !!----------------------------------------------------------------------
 
-      ! Allocate masking array (stored in partition_mod) and workspace array
-      !  for this routine
-      ALLOCATE(imask(jpiglo,jpjglo), zdta(jpiglo,jpjglo), Stat=ierr)
+      ! Allocate masking arrays used in partitioning
+      CALL partition_mask_alloc(jpiglo,jpjglo,ierr)
+      IF(ierr /= 0)THEN
+         CALL ctl_stop('nemo_recursive_partition: failed to allocate masking arrays')
+         RETURN
+      END IF
+
+      ! Allocate local workspace array for this routine
+      ALLOCATE(zdta(jpiglo,jpjglo), Stat=ierr)
       IF(ierr /= 0)THEN
          CALL ctl_stop('nemo_recursive_partition: failed to allocate workspace arrays')
@@ -618 +645 @@
       END IF
 
+      ! Check whether user has specified halo trimming in z via environment variable
+      ! Halo trimming in z is on by default
+      msgtrim_z = .TRUE.
+      CALL GET_ENVIRONMENT_VARIABLE(NAME='NEMO_MSGTRIM_Z', VALUE=lstr, STATUS=ierr)
+      IF( ierr == 0)THEN
+         READ(lstr,FMT="(I)",IOSTAT=ierr) lztrim
+         IF(ierr == 0)THEN
+            IF (lztrim == 0) msgtrim_z = .FALSE.
+         ELSE
+            CALL ctl_warn('nemo_recursive_partition: failed to parse value of NEMO_MSGTRIM_Z environment variable: '//TRIM(lstr))
+         END IF
+      END IF
+
+      WRITE(*,*) 'ARPDBG: msgtrim_z = ',msgtrim_z
+
       ! Factorise the number of MPI PEs to get jpi and jpj as usual
       CALL nemo_partition(num_pes)
 
-      ! Generate a global mask...
-!!$#if defined ARPDEBUG
-!!$      IF(lwp)THEN
-!!$         WRITE(*,*) 'ARPDBG: nemo_recursive_partition: generating mask...'
-!!$         WRITE(*,*) 'ARPDBG: nemo_recursive_partition: jp{i,j}glo = ',jpiglo,jpjglo
-!!$      END IF
-!!$#endif
-
-! ARPDBG - this is the correct variable to check but the dom_nam section
-! of the namelist file hasn't been read in at this stage. 
-!     IF( ntopo == 1 )THEN
-         ! open the file
-         ierr = 0
-!!$         IF ( ln_zco ) THEN 
-!!$            ! Setting ldstop prevents ctl_stop() from being called if the file 
-!!$            ! doesn't exist
-!!$            CALL iom_open ( 'bathy_level.nc', inum, ldstop=.FALSE. ) ! Level bathymetry
-!!$            IF(inum > 0)CALL iom_get ( inum, jpdom_unknown, 'Bathy_level', zdta, &
-!!$                                       kstart=(/jpizoom,jpjzoom/),               &
-!!$                                       kcount=(/jpiglo,jpjglo/) )
-!!$         ELSE
-            CALL iom_open ( 'bathy_meter.nc', inum, ldstop=.FALSE. ) ! Meter bathy in case of partial steps
-            IF(inum > 0)CALL iom_get ( inum, jpdom_unknown, 'Bathymetry' , zdta, &
-                                       kstart=(/jpizoom,jpjzoom/),               &
-                                       kcount=(/jpiglo,jpjglo/) )
-!!$         ENDIF
-         IF(inum > 0)THEN
-            CALL iom_close (inum)
-         ELSE
-            ! Flag that an error occurred when reading the file
-            ierr = 1
-         ENDIF
-!      ELSE
-!         ! Topography not read from file in this case
-!         ierr = 1
-!      END IF
-
-      ! If ln_sco defined then the bathymetry gets smoothed before the 
-      ! simulation begins and that process can alter the coastlines
-      ! therefore we do it here too before calculating our mask.
-!      IF(ln_sco)
-CALL smooth_bathy(zdta)
+      ! ============================
+      ! Generate a global mask from the model bathymetry
+      ! ============================
+
+      ! Read the z-coordinate options from the namelist file
+      REWIND(numnam)
+      READ  (numnam, namzgr)
+
+      ! Read domain options from namelist file
+      CALL dom_nam()
+
+      ! Allocate these arrays so we can use domzgr::zgr_z routine; free them at
+      ! when we're done so as not to upset the 'official' allocation once
+      ! the domain decomposition is done.
+      ALLOCATE(gdepw_0(jpk), gdept_0(jpk), e3w_0(jpk), e3t_0(jpk), &
+               ! Need many global, 3D arrays if zgr_zco is to be called
+               !gdepw(jpiglo,jpjglo,jpk), gdept(jpiglo,jpjglo,jpk), &
+               !gdep3w(jpiglo,jpjglo,jpk), e3t(jpiglo,jpjglo,jpk),  &
+               mig(jpiglo), mjg(jpjglo), &
+               mbathy(jpiglo,jpjglo), bathy(jpiglo,jpjglo), Stat=ierr)
+      IF(ierr /= 0)THEN
+         CALL ctl_stop('nemo_recursive_partition: failed to allocate zgr_z() arrays')
+         RETURN
+      END IF
+
+      ! Set-up reference depth coordinates
+      CALL zgr_z()
+
+      ! Set-up sub-domain limits as global domain for zgr_bat()
+      nldi = 2 ; nlci = jpiglo - 1
+      nldj = 2 ; nlcj = jpjglo - 1
+      jpi = jpiglo
+      jpj = jpjglo
+
+      ! Set-up fake m{i,j}g arrays for zgr_bat() call
+      DO ii = 1, jpiglo, 1
+         mig(ii) = ii
+         mi0(ii) = ii
+         mi1(ii) = ii
+      END DO
+      DO jj = 1, jpjglo, 1
+         mjg(jj) = jj
+         mj0(jj) = jj
+         mj1(jj) = jj
+      END DO
+
+      ! Initialise closed seas so loop over closed seas in zgr_bat works
+      CALL dom_clo()
+
+      ! Read-in bathy (if required) of global domain
+      CALL zgr_bat(.TRUE.)
 
       ! land/sea mask (zero on land, 1 otherwise) over the global/zoom domain
       imask(:,:)=1
-      IF(ierr == 1)THEN
-         ! Failed to read bathymetry so assume all ocean
-         WRITE(*,*) 'ARPDBG: nemo_recursive_partition: no bathymetry file so setting mask to unity'
-
-         ! Mess with otherwise uniform mask to get an irregular decomposition 
-         ! for testing ARPDBG
-         CALL generate_fake_land(imask)
-      ELSE
-         ! Comment-out line below to achieve a regular partition
-         WHERE ( zdta(:,:) <= 1.0E-20 ) imask = LAND
+
+      ! Copy bathymetry in case we need to smooth it
+      zdta(:,:) = bathy(:,:)
+
+      IF(ln_sco)THEN
+         ! If ln_sco defined then the bathymetry gets smoothed before the 
+         ! simulation begins and that process can alter the coastlines (bug!)
+         ! therefore we do it here too before calculating our mask.
+         CALL smooth_global_bathy(zdta, mbathy)
+      ELSE IF(ln_zps)THEN
+         CALL zgr_zps(.TRUE.)
+      ELSE IF(ln_zco)THEN
+         ! Not certain this is required since mbathy computed in zgr_bat()
+         ! in this case.
+         !CALL zgr_zco()
       END IF
+
+      ! Compute the deepest/last ocean level for every point on the grid
+      ibotlevel(:,:) = mbathy(:,:)
+      CALL global_bot_level(ibotlevel)
+
+      ! Comment-out line below to achieve a regular partition
+      WHERE ( zdta(:,:) <= 1.0E-20 ) imask = LAND
 
       ! Allocate partitioning arrays.
@@ -694 +757 @@
 
       ! Now we can do recursive k-section partitioning
-! ARPDBG - BUG if limits on array below are set to anything other than
-! 1 and jp{i,j}glo then check for external boundaries in a few lines
-! time WILL FAIL!
-!      CALL partition_rk ( imask, 1, jpiglo, 1, jpjglo, ierr )
-
-! Multi-core aware version of recursive k-section partitioning
+      ! ARPDBG - BUG if limits on array below are set to anything other than
+      ! 1 and jp{i,j}glo then check for external boundaries in a few lines
+      ! time WILL FAIL!
+      !      CALL partition_rk ( imask, 1, jpiglo, 1, jpjglo, ierr )
+
+      ! Multi-core aware version of recursive k-section partitioning. Currently
+      ! only accounts for whether a grid point is wet or dry. It has no knowledge
+      ! of the number of wet levels at a point.
       CALL partition_mca_rk ( imask, 1, jpiglo, 1, jpjglo, ierr )
 
@@ -708 +773 @@
       ENDIF
 
-      ! Set the mask correctly now we've partitioned
+      ! If we used generate_fake_land() above then we must set
+      ! the mask correctly now we've partitioned. This is only
+      ! necessary when testing.
       !WHERE ( zdta(:,:) <= 0. ) imask = 0
 
-! ARPDBG Quick and dirty dump to stdout in gnuplot form
-!!$      IF(narea == 1)THEN
-!!$         OPEN(UNIT=998, FILE="imask.dat", &
-!!$              STATUS='REPLACE', ACTION='WRITE', IOSTAT=jj)
-!!$         IF( jj == 0 )THEN
-!!$            WRITE (998,*) '# Depth map'
-!!$            DO jj = 1, jpjglo, 1
-!!$               DO ii = 1, jpiglo, 1
-!!$                  WRITE (998,*) ii, jj, zdta(ii,jj) ! imask(ii,jj)
-!!$               END DO
-!!$               WRITE (998,*)
-!!$            END DO
-!!$            CLOSE(998)
-!!$         END IF
-!!$      END IF
+      ! ARPDBG Quick and dirty dump to stdout in gnuplot form
+      IF(narea == 1)THEN
+         OPEN(UNIT=998, FILE="imask.dat", &
+              STATUS='REPLACE', ACTION='WRITE', IOSTAT=jj)
+         IF( jj == 0 )THEN
+            WRITE (998,*) '# Depth map'
+            WRITE (998,*) '# i   j  bathy  imask   ibotlevel   mbathy'
+            DO jj = 1, jpjglo, 1
+               DO ii = 1, jpiglo, 1
+                  WRITE (998,"(I4,1x,I4,1x,E16.6,1x,I4,1x,I4,1x,I4)") &
+                  ii, jj, zdta(ii,jj), imask(ii,jj), ibotlevel(ii,jj), mbathy(ii,jj)
+               END DO
+               WRITE (998,*)
+            END DO
+            CLOSE(998)
+         END IF
+      END IF
 
       jpkm1 = jpk - 1
@@ -742 +811 @@
 
 #if defined ARPDEBUG
+      ! This output is REQUIRED by the check_nemo_comms.pl test script
       WRITE (*,FMT="(I4,' : ARPDBG: ielb, ieub, iesub = ',3I5)") narea-1,&
             ielb, ieub, iesub
@@ -758 +828 @@
       ! false.
       do_exchanges = .TRUE.
+
+      ! Free the domzgr/_oce member arrays that we used earlier in zgr_z() and
+      ! zgr_bat().
+      DEALLOCATE(gdepw_0, gdept_0, e3w_0, e3t_0, mig, mjg,  &
+                 mbathy, bathy)
 
    END SUBROUTINE nemo_recursive_partition

branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/step.F90

@@ -38 +38 @@
    USE asminc           ! assimilation increments    (tra_asm_inc, dyn_asm_inc routines)
    USE timing, ONLY: timing_start, timing_stop, timing_reset, timing_disable
+   USE arpdebugging, ONLY: dump_array
    IMPLICIT NONE
    PRIVATE
@@ -266 +267 @@
       IF(  ln_asmiau .AND. &
          & ln_dyninc       )   CALL dyn_asm_inc( kstp )     ! apply dynamics assimilation increment
+                               !CALL dump_array(kstp,'un_pre_adv',un(:,:,27), &
+                               !                withHalos=.TRUE.)
+
                                !CALL timing_start('dyn_adv') 
                                CALL dyn_adv( kstp )         ! advection (vector or flux form)
                                !CALL timing_stop('dyn_adv','section') 
 
+                               !CALL dump_array(kstp,'ua_pre_vor',ua(:,:,27), &
+                               !                withHalos=.TRUE.)
+
                                !CALL timing_start('dyn_vor') 
                                CALL dyn_vor( kstp )         ! vorticity term including Coriolis
                                !CALL timing_stop('dyn_vor','section') 
 
+                               !CALL dump_array(kstp,'ua_pre_ldf',ua(:,:,27), &
+                               !                withHalos=.TRUE.)
+
                                !CALL timing_start('dyn_ldf') 
                                CALL dyn_ldf( kstp )         ! lateral mixing
@@ -281 +291 @@
       IF(.NOT. Agrif_Root())   CALL Agrif_Sponge_dyn        ! momemtum sponge
 #endif
+                               !CALL dump_array(kstp,'ua_pre_hpg',ua(:,:,27), &
+                               !                withHalos=.TRUE.)
+
                                !CALL timing_start('dyn_hpg') 
                                CALL dyn_hpg( kstp )         ! horizontal gradient of Hydrostatic pressure
                                !CALL timing_stop('dyn_hpg','section') 
 
+                               !CALL dump_array(kstp,'ua_pre_bfr',ua(:,:,27), &
+                               !                withHalos=.TRUE.)
+
                                !CALL timing_start('dyn_bfr') 
                                CALL dyn_bfr( kstp )         ! bottom friction   
                                !CALL timing_stop('dyn_bfr','section') 
+                               
+                               !CALL dump_array(kstp,'ua_pre_zdf',ua(:,:,27), &
+                               !                withHalos=.TRUE.)
 
                                !CALL timing_start('dyn_zdf') 
                                CALL dyn_zdf( kstp )         ! vertical diffusion
                                !CALL timing_stop('dyn_zdf','section') 
+                               
+                               !CALL dump_array(kstp,'ua_pre_spg',ua(:,:,27), &
+                               !                withHalos=.TRUE.)
 
                                !CALL timing_start('dyn_spg') 
                                CALL dyn_spg( kstp, indic )  ! surface pressure gradient
                                !CALL timing_stop('dyn_spg','section') 
+                               !CALL dump_array(kstp,'ua_spg',ua(:,:,27), &
+                               !                withHalos=.TRUE.)
 
                                !CALL timing_start('dyn_nxt') 
                                CALL dyn_nxt( kstp )         ! lateral velocity at next time step
+                               !CALL timing_stop('dyn_nxt','section') 
+
                                !CALL timing_stop('dyn_nxt','section') 
 

branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/timing.F90

@@ -370 +370 @@
       s_wrk => NULL()
       s_timer => s_timer_root
+
+      IF(.NOT. ASSOCIATED(s_timer_root))THEN
+         WRITE(numtime,*) 'No timing information available!'
+         WRITE(numtime,*) '(Have any timed sections been executed?)'
+         RETURN
+      END IF
+
       DO
          ll_ord = .TRUE.
@@ -776 +783 @@
       TYPE(timer), POINTER, INTENT(inout) :: ptr
       !
+      IF(.NOT. ASSOCIATED(ptr))RETURN
+
       IF( ASSOCIATED(ptr%next) ) CALL timing_list(ptr%next)
       IF(lwp) WRITE(numout,*)'   ', ptr%cname   
@@ -816 +825 @@
       !
       TYPE(timer), POINTER :: sl_temp
-    
+ 
+      IF(.NOT. ASSOCIATED(sd_ptr)) RETURN
+
       sl_temp => sd_ptr
       sd_ptr => sd_ptr%next