Index: branches/2011/DEV_r2739_STFC_dCSE/.gitignore
===================================================================
--- branches/2011/DEV_r2739_STFC_dCSE/.gitignore	(revision 3837)
+++ branches/2011/DEV_r2739_STFC_dCSE/.gitignore	(revision 3849)
@@ -6,2 +6,4 @@
 *.err
 *.out
+*.o
+*.exe
Index: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/ARCH/arch-gfortran_linux.fcm
===================================================================
--- branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/ARCH/arch-gfortran_linux.fcm	(revision 3837)
+++ branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/ARCH/arch-gfortran_linux.fcm	(revision 3849)
@@ -15,8 +15,9 @@
 
 
-%NCDF_INC            -I/usr/local/netcdf/include
-%NCDF_LIB            -L/usr/local/netcdf/lib -lnetcdf
+%NCDF_INC            -I/home/kbc59144/MyInstalls/include
+%NCDF_LIB            -L/home/kbc59144/MyInstalls/lib -lnetcdf -lnetcdff
 %FC	             gfortran
-%FCFLAGS             -fdefault-real-8 -O3 -funroll-all-loops -fcray-pointer 
+#%FCFLAGS             -fdefault-real-8 -O3 -funroll-all-loops -fcray-pointer 
+%FCFLAGS             -O0 -g -Wall -Warray-bounds -fbacktrace -fbounds-check -fcheck=all
 %FFLAGS              %FCFLAGS
 %LD                  gfortran
Index: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/ARCH/arch-gfortran_mpi_linux.fcm
===================================================================
--- branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/ARCH/arch-gfortran_mpi_linux.fcm	(revision 3849)
+++ branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/ARCH/arch-gfortran_mpi_linux.fcm	(revision 3849)
@@ -0,0 +1,30 @@
+# generic gfortran compiler options for linux
+# NCDF_INC    netcdf include file
+# NCDF_LIB    netcdf library
+# FC          Fortran compiler command
+# FCFLAGS     Fortran compiler flags
+# FFLAGS      Fortran 77 compiler flags
+# LD          linker
+# LDFLAGS     linker flags, e.g. -L<lib dir> if you have libraries in a
+# FPPFLAGS    pre-processing flags
+# AR          assembler
+# ARFLAGS     assembler flags
+# MK          make
+# USER_INC    additional include files for the compiler,  e.g. -I<include dir>
+# USER_LIB    additional libraries to pass to the linker, e.g. -l<library>
+
+
+%NCDF_INC            -I/home/kbc59144/MyInstalls/include
+%NCDF_LIB            -L/home/kbc59144/MyInstalls/lib -lnetcdf -lnetcdff
+%FC	             mpif90
+%FCFLAGS             -fdefault-real-8 -O3 -funroll-all-loops -fcray-pointer 
+%FFLAGS              %FCFLAGS
+%LD                  mpif90
+%LDFLAGS
+%FPPFLAGS            -P -C -traditional
+%AR                  ar
+%ARFLAGS             -rs
+%MK                  gmake
+%USER_INC            %NCDF_INC
+%USER_LIB            %NCDF_LIB
+
Index: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/CONFIG/ORCA2_LIM/cpp_ORCA2_LIM.fcm
===================================================================
--- branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/CONFIG/ORCA2_LIM/cpp_ORCA2_LIM.fcm	(revision 3837)
+++ branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/CONFIG/ORCA2_LIM/cpp_ORCA2_LIM.fcm	(revision 3849)
@@ -1,1 +1,1 @@
- bld::tool::fppkeys key_trabbl key_vectopt_loop key_orca_r2 key_lim2 key_dynspg_flt key_diaeiv key_ldfslp key_traldf_c2d key_traldf_eiv key_dynldf_c3d key_dtatem key_dtasal key_tradmp key_zdftke key_zdfddm key_zdftmx key_iomput 
+ bld::tool::fppkeys key_trabbl  key_orca_r2 key_lim2 key_dynspg_flt key_diaeiv key_ldfslp key_traldf_c2d key_traldf_eiv key_dynldf_c3d key_dtatem key_dtasal key_tradmp key_zdftke key_zdfddm key_zdftmx   key_mpp_mpi     key_mpp_rkpart
Index: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90
===================================================================
--- branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90	(revision 3837)
+++ branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90	(revision 3849)
@@ -541,7 +541,7 @@
          WHERE( gphit(:,:) < -30._wp)   ;   btm30(:,:) = 0._wp      ! mask out Southern Ocean
 #if defined key_z_first
-         ELSE WHERE                     ;   btm30(:,:) = tmask_1(:,:)
+         ELSEWHERE                     ;   btm30(:,:) = tmask_1(:,:)
 #else
-         ELSE WHERE                     ;   btm30(:,:) = tmask(:,:,1)
+         ELSEWHERE                     ;   btm30(:,:) = tmask(:,:,1)
 #endif
          END WHERE
Index: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/LBC/lib_mpp.F90
===================================================================
--- branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/LBC/lib_mpp.F90	(revision 3837)
+++ branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/LBC/lib_mpp.F90	(revision 3849)
@@ -157,5 +157,6 @@
    INTEGER, SAVE, PUBLIC                 :: nn_cpnode = 4 ! Number of cores per 
                                                           ! compute node on current computer
- 
+   LOGICAL, SAVE, PUBLIC                 :: nn_readpart = .FALSE. ! Whether to read partition from
+                                                            ! file (1) or not (0)
    REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE :: tampon  ! buffer in case of bsend
 
@@ -263,5 +264,6 @@
       !
       NAMELIST/nammpp/ cn_mpi_send, nn_buffer, jpni, jpnj, jpnij, &
-                       nn_xfactors, nn_yfactors, nn_pttrim, nn_cpnode
+                       nn_readpart, nn_xfactors, nn_yfactors,     &
+                       nn_pttrim, nn_cpnode
       !!----------------------------------------------------------------------
       !
Index: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/LBC/mapcomm_mod.F90
===================================================================
--- branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/LBC/mapcomm_mod.F90	(revision 3837)
+++ branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/LBC/mapcomm_mod.F90	(revision 3849)
@@ -156,5 +156,5 @@
 
   ! Public routines
-  PUBLIC :: mapcomms, iprocmap
+  PUBLIC :: mapcomms, iprocmap, set_num_subdomains
 
   ! Public variables
@@ -165,10 +165,9 @@
   PUBLIC :: nsendp,nsendp2d,nrecvp,nrecvp2d,npatchsend,npatchrecv, &
             nxsendp,nysendp,nzsendp,nxrecvp,nyrecvp,nzrecvp,       &
-            idesrecvp,jdesrecvp,isrcsendp,jsrcsendp
+            idesrecvp,jdesrecvp,isrcsendp,jsrcsendp, nprocp
 
   PUBLIC :: ielb,  ieub,  pielb, pjelb, pieub, pjeub,                    &
             iesub, jesub, jeub, ilbext, iubext, jubext, jlbext, pnactive,&
-            piesub, pjesub, jelb, pilbext, pjlbext, pjubext, piubext,    &
-            nprocp
+            piesub, pjesub, jelb, pilbext, pjlbext, pjubext, piubext
 
   PUBLIC :: NONE         &
@@ -191,8 +190,20 @@
 
   ! Switch for trimming points below ocean floor from halo swaps
-  !LOGICAL, PARAMETER :: msgtrim_z = .TRUE. ! .FALSE.
-  LOGICAL, PUBLIC, SAVE      :: msgtrim_z
+  ! Defaults to true unless set via NEMO_MSGTRIM_Z environment var.
+  LOGICAL, PUBLIC, SAVE :: msgtrim_z
 
 CONTAINS
+
+  SUBROUTINE set_num_subdomains(npes)
+    IMPLICIT none
+    !!------------------------------------------------------------------
+    !! Set the number of subdomains to partition the domain into.
+    !!------------------------------------------------------------------
+    INTEGER, INTENT(in) :: npes
+
+    nprocp = npes
+
+  END SUBROUTINE set_num_subdomains
+
 
   SUBROUTINE mapcomms ( depth, ibotlvl, nx, ny, jperio, ierr )
@@ -213,6 +224,6 @@
 
     ! Local variables.
-    INTEGER :: i, i1, i2, icol, ihalo, iproc, iprocc, iprocx, &
-               iprocy, j, j1, j2, lumapout, nadd, naddmaxr, naddmaxs
+    INTEGER :: i, i1, i2, ihalo, iproc, iprocc, iprocx, &
+               iprocy, j, j1, j2, nadd, naddmaxr, naddmaxs
     INTEGER :: ldiff0, ldiff1 ! Local vars for coping with wrapping of coords
     INTEGER :: imax, imin ! Max/min value of i that a halo strip can run 
@@ -1302,5 +1313,5 @@
         IF ( addcorner ) THEN
 #if defined ARPDEBUG
-          WRITE (*,FMT="(I3,': ARPDBG adding corner send to ',I2,', dir = ',I1)") &
+          WRITE (*,FMT="(I3,': ARPDBG adding corner send to ',I4,', dir = ',I1)") &
                  narea-1, procid(iprocc),i
 #endif
@@ -1314,5 +1325,5 @@
 
 #if defined ARPDEBUG
-          WRITE (*,FMT="(I3,': ARPDBG adding corner recv. from ',I3,', old dir = ',I1,' new dir = ',I1)") &
+          WRITE (*,FMT="(I3,': ARPDBG adding corner recv. from ',I4,', old dir = ',I1,' new dir = ',I1)") &
                  narea-1, procid(iprocc),i, j
 #endif
@@ -1380,5 +1391,5 @@
 ! ARP - for debugging only
 !!$        IF(iprocmap == 0)THEN
-!!$           WRITE(*,"('iprocmap: failed to find owner PE for (',I3,I3,')')") ia, ja
+!!$           WRITE(*,"('iprocmap: failed to find owner PE for (',I3,1x,I3,')')") ia, ja
 !!$           WRITE(*,*) 'PE domains are [xmin:xmax][ymin:ymax]:'
 !!$           DO iproc=1,nprocp,1
@@ -1826,5 +1837,5 @@
 
         ! Can any points along the left (low i) edge be trimmed?
-        left_edge: DO i=ilo, ihi - nextra
+        left_edge: DO i=ilo, ihi
           DO j=jlo, jhi
              ! depth is global mask, i and j are local coords
@@ -1856,5 +1867,5 @@
 
         ! Can any points along the right (high i) edge be trimmed?
-        right_edge: DO i=ihi, ilo + nextra, -1
+        right_edge: DO i=ihi, ilo, -1
           DO j=jlo, jhi
 !            IF (depth(i+nimpp-1-nextra,j+njmpp-1) .ne. land) exit right_edge
Index: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/LBC/partition_mod.F90
===================================================================
--- branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/LBC/partition_mod.F90	(revision 3837)
+++ branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/LBC/partition_mod.F90	(revision 3849)
@@ -1,4 +1,4 @@
 MODULE partition_mod
-   USE par_oce, ONLY: jpni, jpnj, jpi, jpj, jpim1, jpjm1, jpij, &
+   USE par_oce, ONLY: jpni, jpnj, jpnij, jpi, jpj, jpim1, jpjm1, jpij, &
                       jpreci, jprecj, jpk, jpkm1, jperio, jpiglo, jpjglo
    USE dom_oce, ONLY: ln_zco, nbondi, nbondj, nidom, npolj, &
@@ -18,5 +18,5 @@
    USE lib_mpp,        ONLY: mppsize, mppsync, mpi_comm_opa,                &
                              MAX_FACTORS, xfactors, yfactors, nn_pttrim,    &
-                             nn_cpnode
+                             nn_cpnode, nn_readpart
 #endif
    USE lib_mpp,        ONLY: ctl_stop, ctl_warn
@@ -34,5 +34,7 @@
                                                  ! (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
+   ! Holds the bottom level of the ocean at each grid point - used for 
+   ! trimming halos in z direction
+   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:,:), TARGET :: ibotlevel 
 
    ! Parameters for the cost function used when evaluating different 
@@ -75,4 +77,5 @@
    PUBLIC imask, ibotlevel, smooth_global_bathy, global_bot_level, partition_mask_alloc
    PUBLIC mpp_init3, partition_rk, partition_mca_rk, write_partition_map
+   PUBLIC read_partition, write_partition
 
 CONTAINS
@@ -292,4 +295,8 @@
       CALL mpp_init_ioipsl()
 
+      ! Write this partition to file in the format that the code can
+      ! read
+      CALL write_partition()
+
       ! ARPDBG - test comms setup
       CALL mpp_test_comms(imask, ibotlevel)
@@ -396,8 +403,4 @@
 #if defined key_mpp_mpi
 
-      ! IMPORTANT: Set the number of PEs to partition over (mapcomm_mod 
-      ! module variable)
-      nprocp = mppsize
-
 #if defined PARTIT_DEBUG
       IF(lwp)WRITE(*,*) 'ARPDBG partition_rk: jpn{i,j} = ',jpni,jpnj
@@ -478,8 +481,4 @@
        ! set that here
        myinst = narea - 1
-
-       ! IMPORTANT: set the number of PEs to partition over (mapcomm_mod 
-       ! module variable)
-       nprocp = mppsize
 
        ! Factorise the total number of MPI processes that we have
@@ -682,4 +681,5 @@
       jpni = nprocx
       jpnj = nprocy
+      jpnij = jpni*jpnj
 
       IF (lwp) THEN
@@ -1388,32 +1388,42 @@
 
 
-   SUBROUTINE finish_partition()
-      USE mapcomm_mod, ONLY: ielb,ieub,pielb,pjelb,pieub,pjeub, &
-                             iesub,jesub,jeub,ilbext,iubext,jubext, &
+   SUBROUTINE finish_partition(fromFile)
+      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, &
+                             piubext, pjlbext, pjubext,                  &
                              trimmed, nidx,eidx,sidx,widx
       IMPLICIT NONE
+      LOGICAL, INTENT(in), OPTIONAL :: fromFile
+      ! Locals
       INTEGER :: iproc, ierr
-
-      ! Set the external boundary flags before boundaries are
-      ! altered by the trimming process and it becomes more difficult
-      ! to recognize which were the external boundaries.
+      LOGICAL :: lFromFile
+
+      ! Check to see whether we're dealing with a partion that has been
+      ! read from file. If we are then there are some things we don't
+      ! calculate here.
+      lFromFile = .FALSE.
+      IF( PRESENT(fromFile) ) lFromFile = fromFile
+
+      IF(.NOT. lFromFile)THEN
+         ! Set the external boundary flags before boundaries are
+         ! altered by the trimming process and it becomes more difficult
+         ! to recognize which were the external boundaries.
       
-      DO iproc=1, nprocp, 1
-           pilbext(iproc) = pielb(iproc).EQ.1
-           piubext(iproc) = pieub(iproc).EQ.jpiglo
-           pjlbext(iproc) = pjelb(iproc).EQ.1
-           pjubext(iproc) = pjeub(iproc).EQ.jpjglo
-        ENDDO
-
-      ! Trim off redundant rows and columns containing all land.
-        IF(.NOT. ALLOCATED(trimmed) )THEN
-           ALLOCATE(trimmed(4,nprocp), Stat=ierr)
-           IF(ierr /= 0)THEN
-              CALL ctl_stop('STOP',    &
-                            'Failed to allocate memory for domain trimming')
-           END IF
-        END IF
+         DO iproc=1, nprocp, 1
+            pilbext(iproc) = pielb(iproc).EQ.1
+            piubext(iproc) = pieub(iproc).EQ.jpiglo
+            pjlbext(iproc) = pjelb(iproc).EQ.1
+            pjubext(iproc) = pjeub(iproc).EQ.jpjglo
+         ENDDO
+
+         ! Trim off redundant rows and columns containing all land.
+         IF(.NOT. ALLOCATED(trimmed) )THEN
+            ALLOCATE(trimmed(4,nprocp), Stat=ierr)
+            IF(ierr /= 0)THEN
+               CALL ctl_stop('STOP',    &
+                             'Failed to allocate memory for domain trimming')
+            END IF
+         END IF
 
 #if defined key_mpp_mpi
@@ -1430,49 +1440,50 @@
         ENDIF
 #else
-        trimmed(1:4,1:nprocp) = .FALSE.
-#endif
-
-        ! Lower boundary (long.) of sub-domain, GLOBAL coords
-        ! before correction for global halos
-        nimpp = pielb(narea)
-
-        ! Is the domain on an external LONGITUDE boundary?
-        nbondi = 0
-        ilbext = pilbext(narea)
-        IF(ilbext)THEN
-           nbondi = -1
-        END IF
-
-        IF( (.NOT. ilbext) .OR. (ilbext .AND. trimmed(widx,narea)) )THEN
-           ! It isn't, which means we must shift its lower boundary by
-           ! -jpreci to allow for the overlap of this domain with its
-           ! westerly neighbour.
-           nimpp = nimpp - jpreci
-        END IF
-
-        iubext = piubext(narea)
-        IF(iubext)THEN
-           nbondi = 1
-           IF(ilbext)nbondi = 2 ! Both East and West boundaries are at 
-                                ! edges of global domain
-        END IF
-
-        ! Set local values for limits in global coords of the sub-domain 
-        ! owned by this process.
-        ielb   = pielb (narea)
-        ieub   = pieub (narea)
-        iesub  = piesub(narea)
-
-        jpi  = iesub + 2*jpreci ! jpi is the same for all domains - this is
-                                ! what original decomposition did
-        nlci = jpi
-
-        ! If the domain is at the edge of the model domain and a cyclic 
-        ! East-West b.c. is in effect then it already incorporates one 
-        ! extra halo column (because of the way the model domain itself is
-        ! set-up). If we've trimmed-off dry rows and columns then, even if
-        ! a domain is on the model boundary, it may still need a halo so
-        ! we add one.
-        IF( nbondi == -1 .AND. (.NOT. trimmed(widx,narea))  )THEN
+         trimmed(1:4,1:nprocp) = .FALSE.
+#endif
+      END IF ! not read from file
+
+      ! Lower boundary (long.) of sub-domain, GLOBAL coords
+      ! before correction for global halos
+      nimpp = pielb(narea)
+
+      ! Is the domain on an external LONGITUDE boundary?
+      nbondi = 0
+      ilbext = pilbext(narea)
+      IF(ilbext)THEN
+         nbondi = -1
+      END IF
+
+      IF( (.NOT. ilbext) .OR. (ilbext .AND. trimmed(widx,narea)) )THEN
+         ! It isn't, which means we must shift its lower boundary by
+         ! -jpreci to allow for the overlap of this domain with its
+         ! westerly neighbour.
+         nimpp = nimpp - jpreci
+      END IF
+
+      iubext = piubext(narea)
+      IF(iubext)THEN
+         nbondi = 1
+         IF(ilbext)nbondi = 2 ! Both East and West boundaries are at 
+                              ! edges of global domain
+      END IF
+
+      ! Set local values for limits in global coords of the sub-domain 
+      ! owned by this process.
+      ielb   = pielb (narea)
+      ieub   = pieub (narea)
+      iesub  = piesub(narea)
+
+      jpi  = iesub + 2*jpreci ! jpi is the same for all domains - this is
+                              ! what original decomposition did
+      nlci = jpi
+
+      ! If the domain is at the edge of the model domain and a cyclic 
+      ! East-West b.c. is in effect then it already incorporates one 
+      ! extra halo column (because of the way the model domain itself is
+      ! set-up). If we've trimmed-off dry rows and columns then, even if
+      ! a domain is on the model boundary, it may still need a halo so
+      ! we add one.
+      IF( nbondi == -1 .AND. (.NOT. trimmed(widx,narea))  )THEN
            ! Western boundary
            ! First column of global domain is actually a halo but NEMO
@@ -1547,8 +1558,8 @@
         nlcj = jpj
 
-! Unlike the East-West boundaries, the global domain does not include
-! halo (rows) at the Northern and Southern edges. In fact, there is no
-! cyclic boundary condition in the North-South direction so there are no
-! halos at all at the edges of the global domain.
+      ! Unlike the East-West boundaries, the global domain does not include
+      ! halo (rows) at the Northern and Southern edges. In fact, there is no
+      ! cyclic boundary condition in the North-South direction so there are no
+      ! halos at all at the edges of the global domain.
       IF( nbondj == -1 .AND. (.NOT. trimmed(sidx,narea)) )THEN
          ! Southern edge
@@ -2772,3 +2783,139 @@
     END SUBROUTINE global_bot_level
 
+
+    SUBROUTINE read_partition(ierr)
+      USE par_oce, ONLY: jpni, jpnj, jpnij
+      USE mapcomm_mod, ONLY: eidx, widx, sidx, nidx, trimmed, &
+                             pilbext, piubext, pjlbext, pjubext
+      IMPLICIT none
+      INTEGER, INTENT(out) :: ierr
+      ! Locals
+      INTEGER, PARAMETER :: funit = 1099
+      INTEGER :: idom, ndom
+      CHARACTER(len=200) :: linein
+      !======================================================================
+
+      ierr = 0
+
+      OPEN(UNIT=funit, file='partition.dat', status='OLD', &
+           ACTION='READ', IOSTAT=ierr)
+      IF(ierr /= 0)THEN
+         CALL ctl_warn('read_partition: failed to read partitioning from file - will calculate it instead.')
+         RETURN
+      END IF
+
+      ! Number of procs in x and y
+      CALL read_next_line(funit, linein, ierr)
+      READ(linein,FMT=*) jpni, jpnj
+
+      ! Store their product
+      jpnij = jpni*jpnj
+
+      ! Check that the implied number of PEs matches that
+      ! in our MPI world
+      ndom = jpni*jpnj
+      IF(ndom /= mppsize)THEN
+         CALL ctl_stop('STOP', &
+                       'read_partition: no. of PEs specified in partition.dat does not match no. of PEs in use by this job.')
+      END IF
+
+      ! Read the description of each sub-domain
+      domains: DO idom = 1, ndom, 1
+
+         ! Coordinates of bottom-left (SW) corner of domain
+         CALL read_next_line(funit, linein, ierr)
+         READ(linein,FMT=*) pielb(idom), pjelb(idom)
+         ! Top-right (NE) corner
+         CALL read_next_line(funit, linein, ierr)
+         READ(linein,FMT=*) pieub(idom), pjeub(idom)
+         ! Whether this domain has external boundaries and has been trimmed
+         CALL read_next_line(funit, linein, ierr)
+         READ(linein,FMT=*) pilbext(idom), trimmed(widx,idom)
+         CALL read_next_line(funit, linein, ierr)
+         READ(linein,FMT=*) piubext(idom), trimmed(eidx,idom)
+         CALL read_next_line(funit, linein, ierr)
+         READ(linein,FMT=*) pjlbext(idom), trimmed(sidx,idom)
+         CALL read_next_line(funit, linein, ierr)
+         READ(linein,FMT=*) pjubext(idom), trimmed(nidx,idom)
+
+         piesub(idom) = pieub(idom) - pielb(idom) + 1
+         pjesub(idom) = pjeub(idom) - pjelb(idom) + 1
+
+      END DO domains
+
+      ! All done - close the file
+      CLOSE(UNIT=funit)
+
+      CALL finish_partition(fromFile=.TRUE.)
+
+    END SUBROUTINE read_partition
+
+    !===================================================================
+
+    SUBROUTINE write_partition
+      USE par_oce, ONLY: jpni, jpnj
+      USE mapcomm_mod, ONLY: eidx, widx, sidx, nidx, trimmed,    &
+                             pjubext, pjlbext, piubext, pilbext, &
+                             pielb, pieub, pjelb, pjeub
+      IMPLICIT none
+      INTEGER, PARAMETER :: funit = 1099
+      INTEGER :: ierr
+      INTEGER :: idom
+
+      ! Only PE 0 (narea==1) writes this file
+      IF(narea /= 1) RETURN
+
+      OPEN(UNIT=funit, file='partition.dat.new', status='REPLACE', &
+           ACTION='WRITE', IOSTAT=ierr)
+      IF(ierr /= 0)THEN
+         CALL ctl_warn('write_partition: failed to write partition description to file.')
+         RETURN
+      END IF
+      WRITE(UNIT=funit,FMT="('#  jpni  jpnj')")
+      WRITE(UNIT=funit,FMT="(I5,1x,I5)") jpni, jpnj
+
+      DO idom = 1, mppsize, 1
+         WRITE(UNIT=funit,FMT="('# Domain: ',I5)") idom
+         IF(idom==1)WRITE(UNIT=funit,FMT="('# Lower bounds: x  y')")
+         WRITE(UNIT=funit,FMT="(I5,1x,I5)") pielb(idom), pjelb(idom)
+         IF(idom==1)WRITE(UNIT=funit,FMT="('# Upper bounds: x  y')")
+         WRITE(UNIT=funit,FMT="(I5,1x,I5)") pieub(idom), pjeub(idom)
+         IF(idom==1)WRITE(UNIT=funit,FMT="('# x: Lower bound external, trimmed')")
+         WRITE(UNIT=funit,FMT="(L5,1x,L5)") pilbext(idom), trimmed(widx,idom)
+         IF(idom==1)WRITE(UNIT=funit,FMT="('# x: Upper bound external, trimmed')")
+         WRITE(UNIT=funit,FMT="(L5,1x,L5)") piubext(idom), trimmed(eidx,idom)
+         IF(idom==1)WRITE(UNIT=funit,FMT="('# y: Lower bound external, trimmed')")
+         WRITE(UNIT=funit,FMT="(L5,1x,L5)") pjlbext(idom), trimmed(sidx,idom)
+         IF(idom==1)WRITE(UNIT=funit,FMT="('# y: Upper bound external, trimmed')")
+         WRITE(UNIT=funit,FMT="(L5,1x,L5)") pjubext(idom), trimmed(nidx,idom)
+      END DO
+
+      CLOSE(UNIT=funit)
+
+    END SUBROUTINE write_partition
+
+    SUBROUTINE read_next_line(funit, linein, ierr)
+      IMPLICIT none
+      !!------------------------------------------------------------------
+      INTEGER,            INTENT( in) :: funit  ! Unit no. to read
+      CHARACTER(len=200), INTENT(out) :: linein ! String containing next
+                                                ! non-comment line in file
+      INTEGER,            INTENT(out) :: ierr   ! Error flag (0==OK)
+      !!------------------------------------------------------------------
+
+      ierr = 0
+
+      READ(UNIT=funit,FMT="(200A)") linein
+
+      ! Comment lines begin with '#'. Skip those plus any blank
+      ! lines...
+      DO WHILE( INDEX( TRIM(ADJUSTL(linein)),'#') /= 0 .OR. &
+                LEN_TRIM(linein) == 0 )
+         READ(UNIT=funit,FMT="(200A)") linein   
+      END DO
+
+      WRITE(*,*)'returning linein >>'//linein//'<<'
+
+    END SUBROUTINE read_next_line
+
 END MODULE partition_mod
Index: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/OBS/obs_readmdt.F90
===================================================================
--- branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/OBS/obs_readmdt.F90	(revision 3837)
+++ branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/OBS/obs_readmdt.F90	(revision 3849)
@@ -114,5 +114,5 @@
       WHERE(z_mdt(:,:) /= zfill)   ;   mdtmask(:,:) = tmask(:,:,1)
 #endif
-      ELSE WHERE                   ;   mdtmask(:,:) = 0
+      ELSEWHERE                   ;   mdtmask(:,:) = 0
       END WHERE
 
Index: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90
===================================================================
--- branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90	(revision 3837)
+++ branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90	(revision 3849)
@@ -105,13 +105,13 @@
       ! ice free of snow and melts
       WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalbfz(:,:,:) = rn_albice
-      ELSE WHERE                                              ;   zalbfz(:,:,:) = rn_alphdi
+      ELSEWHERE                                              ;   zalbfz(:,:,:) = rn_alphdi
       END  WHERE
 
       WHERE     ( 1.5  < ph_ice                     )  ;  zficeth = zalbfz
-      ELSE WHERE( 1.0  < ph_ice .AND. ph_ice <= 1.5 )  ;  zficeth = 0.472  + 2.0 * ( zalbfz - 0.472 ) * ( ph_ice - 1.0 )
-      ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.0 )  ;  zficeth = 0.2467 + 0.7049 * ph_ice              &
+      ELSEWHERE( 1.0  < ph_ice .AND. ph_ice <= 1.5 )  ;  zficeth = 0.472  + 2.0 * ( zalbfz - 0.472 ) * ( ph_ice - 1.0 )
+      ELSEWHERE( 0.05 < ph_ice .AND. ph_ice <= 1.0 )  ;  zficeth = 0.2467 + 0.7049 * ph_ice              &
          &                                                                 - 0.8608 * ph_ice * ph_ice     &
          &                                                                 + 0.3812 * ph_ice * ph_ice * ph_ice
-      ELSE WHERE                                       ;  zficeth = 0.1    + 3.6    * ph_ice
+      ELSEWHERE                                       ;  zficeth = 0.1    + 3.6    * ph_ice
       END WHERE
 
Index: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90
===================================================================
--- branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90	(revision 3837)
+++ branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90	(revision 3849)
@@ -1433,5 +1433,5 @@
       !
       INTEGER ::   jit, jk   ! dummy loop indices
-      INTEGER ::   id1, id2, id3, id4, id5, id6, id7, id8
+      INTEGER ::   id1, id2, id3, id4, id5, id6
       INTEGER ::   ji, jj, ikbu, ikbv
       REAL(wp)::   cbx, cby
@@ -1448,5 +1448,5 @@
             id6 = iom_varid( numror, 'mxln' , ldstop = .FALSE. )
             !
-            IF( MIN( id1, id2, id3, id4, id5, id6, id7, id8 ) > 0 ) THEN        ! all required arrays exist
+            IF( MIN( id1, id2, id3, id4, id5, id6 ) > 0 ) THEN        ! all required arrays exist
                CALL iom_get( numror, jpdom_autoglo, 'en'    , en     )
                CALL iom_get( numror, jpdom_autoglo, 'avt'   , avt    )
Index: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
===================================================================
--- branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 3837)
+++ branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 3849)
@@ -601,10 +601,11 @@
       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, &
+      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, msgtrim_z
-      USE partition_mod, ONLY: partition_rk, partition_mca_rk, &
-                               imask, ibotlevel, partition_mask_alloc, &
+                             piubext, pjlbext, pjubext, LAND, trimmed,        &
+                             msgtrim_z, set_num_subdomains
+      USE partition_mod, ONLY: partition_rk, partition_mca_rk, read_partition, &
+                               imask, ibotlevel, partition_mask_alloc,         &
                                smooth_global_bathy, global_bot_level
       USE par_oce,       ONLY: do_exchanges
@@ -645,8 +646,10 @@
       END IF
 
-      ! Check whether user has specified halo trimming in z via environment variable
+      ! 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)
+      CALL GET_ENVIRONMENT_VARIABLE(NAME='NEMO_MSGTRIM_Z', VALUE=lstr, &
+                                    STATUS=ierr)
       IF( ierr == 0)THEN
          READ(lstr,FMT="(I)",IOSTAT=ierr) lztrim
@@ -658,8 +661,5 @@
       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)
+      IF(lwp) WRITE(*,*) 'ARPDBG: msgtrim_z = ',msgtrim_z
 
       ! ============================
@@ -674,15 +674,13 @@
       CALL dom_nam()
 
-      ! Allocate these arrays so we can use domzgr::zgr_z routine; free them at
+      ! Allocate these arrays so we can use domzgr::zgr_z routine; free them
       ! 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')
+         CALL ctl_stop('STOP', &
+                       'nemo_recursive_partition: failed to allocate zgr_z() arrays')
          RETURN
       END IF
@@ -743,10 +741,10 @@
       ! Allocate partitioning arrays.
 
-      IF ( .not.allocated(pielb) ) THEN
+      IF ( .NOT. ALLOCATED(pielb) ) THEN
          ALLOCATE (pielb(num_pes),   pieub(num_pes), piesub(num_pes),     &
                    pilbext(num_pes), piubext(num_pes),                    &
                    pjelb(num_pes),   pjeub(num_pes), pjesub(num_pes),     &
                    pjlbext(num_pes), pjubext(num_pes), pnactive(num_pes), &
-                   Stat = ierr)
+                   trimmed(4,num_pes), Stat = ierr)
          IF(ierr /= 0)THEN
             CALL ctl_stop('STOP', &
@@ -756,14 +754,32 @@
       ENDIF
 
-      ! 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. 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 )
+      ! Set error flag so that we calculate domain decomp if not reading
+      ! existing decomposition or if read fails.
+      ierr = 1
+
+      IF( nn_readpart )THEN
+         ! Read the partitioning to use from disk 
+         CALL read_partition(ierr)
+         IF ( ierr /= 0 ) THEN
+            CALL ctl_warn('Read of pre-calculated domain decomposition failed - will calculate one instead.')
+         END IF
+      END IF
+
+      ! Set the number of sub-domains for which we are to partition 
+      ! (module var in mapcomm_mod)
+      CALL set_num_subdomains(num_pes)
+
+      IF(ierr /= 0)THEN
+         ! 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 )
+
+         ! 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 )
+     END IF
 
       ! Check the error code from partitioning.
@@ -854,5 +870,5 @@
       ! Search backwards from the square root of n.
 
-      fact_loop: DO kna=SQRT(REAL(kn)),1,-1
+      fact_loop: DO kna=INT(SQRT(REAL(kn))),1,-1
          IF ( kn/kna*kna == kn ) THEN
             EXIT fact_loop
Index: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/TOOLS/RK_PARTITION/src/make_partition.f90
===================================================================
--- branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/TOOLS/RK_PARTITION/src/make_partition.f90	(revision 3849)
+++ branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/TOOLS/RK_PARTITION/src/make_partition.f90	(revision 3849)
@@ -0,0 +1,173 @@
+PROGRAM make_partition
+   USE iso_c_binding ! For command-line interface
+   USE partition_mod, ONLY: jpiglo, jpjglo, write_partition, nemo_partition, &
+                            allocate_arrays, partition_mca_rk, cyclic_bc
+   USE mask_mod, ONLY: imask, generate_mask
+   IMPLICIT none
+
+   INTEGER :: narg, iarg, arglen
+   INTEGER :: ji, ierr
+   INTEGER :: num_pes
+   CHARACTER(len=100) :: argval
+   INTEGER :: cores_per_node ! No. of cores per compute node on the machine
+                             ! that NEMO will run on.
+   LOGICAL :: regular ! True for regular decomp, false otherwise
+
+   !!----------------------------------------------------------------------
+
+   ! Initialise domain size
+   jpiglo  = 0 ; jpjglo = 0
+   num_pes = 0
+   cores_per_node = 16
+   regular   = .FALSE.        ! By default we do recursive partitioning
+   cyclic_bc = .FALSE.        ! Cyclic boundary conditions off by default
+
+   ! Parse the command-line arguments
+   narg = COMMAND_ARGUMENT_COUNT()
+
+   IF(narg == 0)THEN
+      CALL print_usage()
+      STOP
+   END IF
+
+   iarg = 1
+   parse_args: DO WHILE(iarg <= narg)
+
+      CALL GET_COMMAND_ARGUMENT(iarg,argval,arglen,ierr)
+
+      IF(TRIM(argval) == '-d')THEN
+         iarg = iarg + 1
+         CALL GET_COMMAND_ARGUMENT(iarg,argval,arglen,ierr)
+
+         ! Parse this argument. We expect it to be of the form <nx>x<ny>.
+         ji = INDEX(argval,"x")
+         IF(ji < 2 .OR. ji > (arglen-1))THEN
+            WRITE(*,"('ERROR: unrecognised argument format: ',(A))") TRIM(argval)
+         END IF
+         READ(argval(1:ji-1),*,IOSTAT=ierr) jpiglo
+         IF(ierr /= 0)THEN
+            WRITE(*,"('ERROR: failed to parse jpiglo argument')")
+            STOP
+         END IF
+         READ(argval(ji+1:),*,IOSTAT=ierr) jpjglo
+         IF(ierr /= 0)THEN
+            WRITE(*,"('ERROR: failed to parse jpjglo argument')")
+            STOP
+         END IF
+
+      ELSE IF(TRIM(argval) == '-np')THEN
+         iarg = iarg + 1
+         CALL GET_COMMAND_ARGUMENT(iarg,argval,arglen,ierr)
+
+         ! Parse this argument. We expect it to be of the form <np>.
+         READ(argval,*,IOSTAT=ierr) num_pes
+         IF(ierr /= 0)THEN
+            WRITE(*,"('ERROR: failed to parse num_pes argument')")
+            STOP
+         END IF
+
+      ELSE IF(TRIM(argval) == '-c')THEN
+         iarg = iarg + 1
+         CALL GET_COMMAND_ARGUMENT(iarg,argval,arglen,ierr)
+         ! Parse this argument. We expect it to be of the form <cpn>.
+         READ(argval,*,IOSTAT=ierr) cores_per_node
+
+      ELSE IF(TRIM(argval) == '-r')THEN
+         ! Switch-on regular partitioning
+         regular = .TRUE.
+
+      ELSE IF(TRIM(argval) == '-cyclic')THEN
+         ! Switch-on cyclic boundary conditions
+         cyclic_bc = .TRUE.
+
+      ELSE IF(TRIM(argval) == '-h')THEN
+         CALL print_usage()
+
+      ELSE
+         WRITE(*,"('ERROR: unrecognised argument: ',(A))") TRIM(argval)
+         STOP
+
+      END IF
+
+      iarg = iarg + 1
+   END DO parse_args
+
+   ! Check that essential parameters have been supplied by user
+   IF( num_pes == 0 )THEN
+      WRITE (*,*) 'ERROR: you must specify num. PEs to decompose over!'
+      CALL print_usage()
+      STOP
+   END IF
+
+   ! Check whether domain dimensions have been specified
+   IF( jpiglo < 1 .OR. jpjglo < 1 )THEN
+      IF( regular )THEN
+         ! Actually, I could read the bathy file if available to get the domain
+         ! dimensions for this case.
+         WRITE (*,*) 'ERROR: you must specify domain extent when doing regular partition!'
+         CALL print_usage()
+         STOP
+      ELSE
+         WRITE(*,*) 'Domain dimensions will be read from bathymetry file'
+      END IF
+   END IF
+
+   WRITE (*,FMT="(/'No. of MPI Processes = ',I4)") num_pes
+   IF(.NOT. regular)THEN
+      WRITE (*,FMT="('No. of cores per compute node = ',I4)") cores_per_node
+   END IF
+
+   ! Generate the land/sea mask from the bathymetry which is read
+   ! from file
+   CALL generate_mask(regular, ierr)
+   IF(ierr /= 0)THEN
+      STOP 'Failed to generate land/sea mask'
+   END IF
+
+   CALL allocate_arrays(num_pes, ierr)
+
+   IF(regular)THEN
+
+      ! Imitate standard NEMO regular partitioning
+      CALL nemo_partition(num_pes)
+
+   ELSE
+
+      ! 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, num_pes, cores_per_node, ierr )
+
+      ! Check the error code from partitioning.
+      IF ( ierr /= 0 ) THEN
+         STOP 'nemo_recursive_partition: Partitioning failed'
+      ENDIF
+
+   END IF
+
+   CALL write_partition(imask)
+
+   WRITE(*,FMT="(/'Partitioning done and written to partition.dat.new. Also'/ &
+        &  'written in postscript to domain_decomp.ps.'/)")
+
+CONTAINS
+
+  SUBROUTINE print_usage()
+     IMPLICIT None
+     WRITE(*,*)
+     WRITE(*,*) 'Usage: make_partition.exe -d <jpiglo>x<jpjglo> -np <nproc> -c <cores_per_node> [-r] [-cyclic] [-h]'
+     WRITE(*,*) 'Where:  domain is jpiglo x jpjglo in extent'
+     WRITE(*,*) '        nproc is no. of sub-domains/MPI procs required'
+     WRITE(*,*) '        cores_per_node is no. of cores per compute node on machine'
+     WRITE(*,*) '           on which NEMO is to be executed'
+     WRITE(*,*) '        -r      : generate regular (old style) partition'
+     WRITE(*,*) '        -cyclic : enable cyclic E/W boundary conditions'
+     WRITE(*,*) '        -h       : print this message'
+     WRITE(*,*)
+     WRITE(*,*) 'Vertical coordinate options are read from the NEMO namelist file and'
+     WRITE(*,*) 'used to determine whether to smooth the bathy.'
+     WRITE(*,*)
+  END SUBROUTINE print_usage
+
+
+END PROGRAM make_partition
Index: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/TOOLS/RK_PARTITION/src/mask_mod.f90
===================================================================
--- branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/TOOLS/RK_PARTITION/src/mask_mod.f90	(revision 3849)
+++ branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/TOOLS/RK_PARTITION/src/mask_mod.f90	(revision 3849)
@@ -0,0 +1,311 @@
+MODULE mask_mod
+  USE par_kind, ONLY: wp
+  IMPLICIT none
+  PRIVATE
+
+  ! To hold the wet/dry mask for the domain. Obtained from the bathymetry
+  ! read from file.
+  INTEGER, ALLOCATABLE, DIMENSION(:,:) :: imask
+
+  LOGICAL :: ln_zco, ln_zps, ln_sco        ! Vert. coordinate options
+
+   !* Namelist namzgr_sco *
+   REAL(wp) ::  rn_sbot_min = 300._wp    ! minimum depth of s-bottom surface (>0) (m)
+   REAL(wp) ::  rn_sbot_max = 5250._wp    ! maximum depth of s-bottom surface (= ocean depth) (>0) (m)
+   REAL(wp) ::  rn_theta    =    6.00_wp  ! surface control parameter (0<=rn_theta<=20)
+   REAL(wp) ::  rn_thetb    =    0.75_wp  ! bottom control parameter  (0<=rn_thetb<= 1)
+   REAL(wp) ::  rn_rmax     =    0.15_wp  ! maximum cut-off r-value allowed (0<rn_rmax<1)
+   LOGICAL  ::  ln_s_sigma  = .false.     ! use hybrid s-sigma -coordinate & stretching function fssig1 (ln_sco=T)
+   REAL(wp) ::  rn_bb       =    0.80_wp  ! stretching parameter for song and haidvogel stretching
+   !                                       ! ( rn_bb=0; top only, rn_bb =1; top and bottom)
+   REAL(wp) ::  rn_hc       =  150._wp    ! Critical depth for s-sigma coordinates
+
+  PUBLIC imask
+  PUBLIC generate_mask
+
+CONTAINS
+
+  SUBROUTINE read_namelist_options(ierr)
+    !!---------------------------------------------------------------------
+    ! Read the z-coordinate options from the namelist file
+    !!---------------------------------------------------------------------
+    IMPLICIT none
+    INTEGER, INTENT(out) :: ierr
+    ! Locals
+    INTEGER, PARAMETER :: numnam = 31
+    ! From OPA_SRC/DOM/domzgr
+    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
+
+    OPEN( UNIT=numnam, FILE='namelist', STATUS='OLD', FORM='FORMATTED', &
+         ACCESS='SEQUENTIAL',IOSTAT=ierr )
+    IF(ierr /= 0)THEN
+       WRITE(*,*) 'read_namelist_options: ERROR: failed to open file ''namelist'''
+    END IF
+      
+    REWIND(numnam)
+    READ  (numnam, namzgr)
+
+    IF(ln_sco)THEN
+       REWIND(numnam)
+       READ(numnam, namzgr_sco)
+    END IF
+
+    CLOSE(numnam)
+
+  END SUBROUTINE read_namelist_options
+
+
+  SUBROUTINE generate_mask(regular, ierr)
+    !!---------------------------------------------------------------------
+    !!---------------------------------------------------------------------
+    USE netcdf
+    USE partition_mod, ONLY: LAND, jpiglo, jpjglo
+    IMPLICIT none
+    LOGICAL, INTENT(In) :: regular
+    INTEGER, INTENT(out) :: ierr
+    ! Local vars
+    INTEGER :: if90id           ! nf90 identifier of the opened file
+    INTEGER :: ivarid           ! netcdf variable ID
+    INTEGER :: i_dimid_lat, i_dimid_lon ! IDs of dimensions in netcdf file
+    INTEGER :: ilat, ilon               ! Current values of those dimensions
+    REAL(KIND(1.0D0)), ALLOCATABLE, DIMENSION(:,:) :: bathy
+    !!---------------------------------------------------------------------  
+
+    ! Initialise error flag
+    ierr = 0
+
+    ! We read the NEMO namelist file to establish what z-coordinate
+    ! options are being used. If ln_sco is .TRUE. then we read
+    ! the parameters used in smoothing the bathymetry.
+    CALL read_namelist_options(ierr)
+    IF( ierr /= 0)THEN
+       WRITE(*,*) 'Failed to read z-coord options from namelist file'
+       RETURN
+    END IF
+
+    IF(.NOT. regular)THEN
+       ! ============================
+       ! Generate a global mask from the model bathymetry
+       ! ============================
+
+       ! Read the bathymetry from file
+       ! Open the file
+       ierr = NF90_OPEN( 'bathy_meter.nc', NF90_NOWRITE  , if90id )
+       IF(ierr /= nf90_noerr)THEN
+          WRITE(*,*) 'generate_mask: Failed to open netcdf file bathy_meter.nc'
+          ierr = 1
+          RETURN
+       END IF
+
+       IF(jpiglo < 1 .OR. jpjglo < 1)THEN
+
+          ierr = NF90_INQ_DIMID(if90id, 'lat', i_dimid_lat )
+          IF(ierr /= nf90_noerr)THEN
+             WRITE(*,*) 'generate_mask: Failed to read ID of lat dimension from netcdf file bathy_meter.nc'
+             ierr = 1
+             RETURN
+          END IF
+
+          ierr = NF90_INQ_DIMID(if90id, 'lon', i_dimid_lon )
+          IF(ierr /= nf90_noerr)THEN
+             WRITE(*,*) 'generate_mask: Failed to read ID of lon dimension from netcdf file bathy_meter.nc'
+             ierr = 1
+             RETURN
+          END IF
+
+          ierr = NF90_INQUIRE_DIMENSION( if90id, i_dimid_lat, len = ilat )
+          IF(ierr /= nf90_noerr)THEN
+             WRITE(*,*) 'generate_mask: Failed to read lat dimension from netcdf file bathy_meter.nc'
+             ierr = 1
+             RETURN
+          END IF
+       
+          ierr = NF90_INQUIRE_DIMENSION( if90id, i_dimid_lon, len = ilon )
+          IF(ierr /= nf90_noerr)THEN
+             WRITE(*,*) 'generate_mask: Failed to read lon dimension from netcdf file bathy_meter.nc'
+             ierr = 1
+             RETURN
+          END IF
+
+          jpiglo = ilon
+          jpjglo = ilat
+
+       END IF ! jpiglo or jpjglo unset
+
+       ! Find the variable ID
+       ierr = NF90_INQ_VARID( if90id, 'Bathymetry', ivarid )
+       IF( ierr /= nf90_noerr)THEN
+          WRITE (*,*) 'generate_mask: ERROR: netcdf file does not contain ''Bathymetry'' variable'
+          ierr = 1
+          RETURN
+       END IF
+
+       WRITE (*,FMT="('Dimensions of domain are ',I4,' x ',I4)") jpiglo,jpjglo
+
+       ! Allocate masking arrays used in partitioning
+       ALLOCATE(bathy(jpiglo,jpjglo), imask(jpiglo,jpjglo), STAT=ierr)
+       IF(ierr /= 0)THEN
+          WRITE (*,*) 'ERROR: failed to allocate memory for bathy and mask arrays'
+          RETURN
+       END IF
+
+       ! Read the whole field - start is 1,1 by default
+       ierr = NF90_GET_VAR(if90id, ivarid, bathy, count=(/jpiglo,jpjglo/) )
+       IF( ierr /= nf90_noerr )THEN
+          WRITE (*,*) 'generate_mask: ERROR: failed to read ''Bathymetry'' variable from netcdf file'
+          ierr = 1
+          RETURN
+       END IF
+
+       ! Smooth the bathymetry if required
+       IF(ln_sco)THEN
+          CALL smooth_bathy(bathy,ierr)
+          IF(ierr /= 0)THEN
+             WRITE(*,*) 'generate_mask: ERROR: failed to smooth bathy'
+             RETURN
+          END IF
+       END IF
+
+       ! Set imask using the bathymetry (zero on land, 1 otherwise) over 
+       ! the global/zoom domain.
+       WHERE ( bathy(:,:) <= 1.0E-20 ) 
+          imask = LAND
+       ELSEWHERE
+          imask = 1
+       END WHERE
+
+       WRITE (*,FMT="(/'generate_mask: Land/sea mask set from bathymetry OK.')")
+
+    END IF
+
+  END SUBROUTINE generate_mask
+
+
+  SUBROUTINE smooth_bathy(inbathy, ierr)
+    IMPLICIT none
+    !!----------------------------------------------------------------------
+    !!                      Routine smooth_bathy
+    !!   Replicates the smoothing done on the decomposed domain in zgr_sco()
+    !!   in domzgr.F90. However, here the domain is NOT decomposed and
+    !!   every processor performs an identical smoothing on the whole model
+    !!   bathymetry. This is to ensure that the domain decomposition
+    !!   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 (erroneously) change whether grid points are wet or dry.)
+    !!----------------------------------------------------------------------
+    REAL(wp), INTENT(inout), DIMENSION(:,:) :: inbathy ! The bathymetry to 
+                                                       ! be smoothed
+    INTEGER, INTENT(out)                    :: ierr
+    ! Locals
+    INTEGER  :: ji, jj, jl
+    INTEGER  :: iip1, ijp1, iim1, ijm1   ! temporary integers
+    INTEGER  :: x_size, y_size
+    REAL(wp) :: zrmax, zri, zrj
+    REAL(wp), PARAMETER :: TOL_ZERO = 1.0E-20_wp ! Any value less than 
+                                                 ! this assumed zero
+    REAL(wp), ALLOCATABLE, DIMENSION(:,:)   :: zenv, ztmp, zmsk, zbot, &
+                                               zscosrf, zhbatt
+    !
+    !!----------------------------------------------------------------------
+
+    ! Do this because we've not decomposed the domain yet and therefore
+    ! jpi,jpj,nlc{i,j} etc. are not set.
+    x_size = SIZE(inbathy, 1)
+    y_size = SIZE(inbathy, 2)
+
+    ALLOCATE(zenv(x_size,y_size), ztmp(x_size,y_size), zmsk(x_size,y_size), &
+             zbot(x_size,y_size), zhbatt(x_size, y_size), &
+             zscosrf(x_size,y_size), Stat=ierr)
+    IF( ierr /= 0 ) THEN
+       WRITE(*,*) 'smooth_bathy: ERROR - failed to allocate workspace arrays'
+       RETURN
+    END IF
+
+    zscosrf(:,:) = 0._wp      ! ocean surface depth (here zero: no under ice-shelf sea)
+    zbot(:,:) = inbathy(:,:)  ! ocean bottom depth
+    !                         ! set maximum ocean depth
+    inbathy(:,:) = MIN( rn_sbot_max, inbathy(:,:) )
+
+    WHERE( inbathy(:,:) > TOL_ZERO ) inbathy(:,:) = &
+                                       MAX( rn_sbot_min, inbathy(:,:) )
+
+    ! use r-value to create hybrid coordinates
+    zenv(:,:) = MAX( inbathy(:,:), rn_sbot_min )
+    ! 
+    ! Smooth the bathymetry (if required)
+    !
+    jl = 0
+    zrmax = 1._wp
+  
+    WRITE(*,"(/'Smoothing bathymetry....'/)")
+                                                       ! ================ !
+    DO WHILE( jl <= 10000 .AND. zrmax > rn_rmax )      !  Iterative loop  !
+       !                                               ! ================ !
+       jl = jl + 1
+       zrmax = 0._wp
+       zmsk(:,:) = 0._wp
+
+       DO jj = 1, y_size
+          DO ji = 1, x_size
+             iip1 = MIN( ji+1, x_size )      ! force zri = 0 on last line (ji=ncli+1 to jpi)
+             ijp1 = MIN( jj+1, y_size )      ! force zrj = 0 on last row  (jj=nclj+1 to jpj)
+             zri = ABS( zenv(iip1,jj  ) - zenv(ji,jj) ) / ( zenv(iip1,jj  ) + zenv(ji,jj) )
+             zrj = ABS( zenv(ji  ,ijp1) - zenv(ji,jj) ) / ( zenv(ji  ,ijp1) + zenv(ji,jj) )
+             zrmax = MAX( zrmax, zri, zrj )
+             
+             IF( zri > rn_rmax )   zmsk(ji  ,jj  ) = 1._wp
+             IF( zri > rn_rmax )   zmsk(iip1,jj  ) = 1._wp
+             IF( zrj > rn_rmax )   zmsk(ji  ,jj  ) = 1._wp
+             IF( zrj > rn_rmax )   zmsk(ji  ,ijp1) = 1._wp
+          END DO
+       END DO
+
+       !
+       WRITE(*,"('smooth_bathy : iter=',I5,' rmax=',F8.4,' nb of pt= ',I8)") &
+                                             jl, zrmax, INT( SUM(zmsk(:,:) ) )
+       !
+
+       ! Copy current surface before next smoothing iteration 
+       ztmp(:,:) = zenv(:,:)
+
+       DO jj = 1, y_size
+          DO ji = 1, x_size
+             iip1 = MIN( ji+1, x_size )     ! last  line (ji=nlci)
+             ijp1 = MIN( jj+1, y_size )     ! last  raw  (jj=nlcj)
+             iim1 = MAX( ji-1,  1  )      ! first line (ji=nlci)
+             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)  &
+             &                      )
+             ENDIF
+          END DO
+       END DO
+       !
+       DO jj = 1,y_size
+          DO ji = 1,x_size
+             IF( zmsk(ji,jj) >= 1._wp-TOL_ZERO ) zenv(ji,jj) = MAX( ztmp(ji,jj), inbathy(ji,jj) )
+          END DO
+       END DO
+       !
+       !                                                ! ================ !
+    END DO                                              !     End loop     !
+    !                                                   ! ================ !
+    !
+    ! envelop bathymetry saved in zhbatt
+    zhbatt(:,:) = zenv(:,:) 
+
+    ! Subtract off rn_sbot_min so can check for land using zenv = LAND (0)
+    inbathy(:,:) = zenv(:,:) - rn_sbot_min
+
+  END SUBROUTINE smooth_bathy
+
+END MODULE mask_mod
Index: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/TOOLS/RK_PARTITION/src/par_kind.F90
===================================================================
--- branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/TOOLS/RK_PARTITION/src/par_kind.F90	(revision 3849)
+++ branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/TOOLS/RK_PARTITION/src/par_kind.F90	(revision 3849)
@@ -0,0 +1,40 @@
+MODULE par_kind
+   !!======================================================================
+   !!                   ***  MODULE par_kind  ***
+   !! Ocean :  define the kind of real for the whole model
+   !!======================================================================
+   !! History :  1.0  ! 2002-06  (G. Madec)  Original code
+   !!            3.3  ! 2010-12  (G. Madec)  add a standard length of character strings
+   !!----------------------------------------------------------------------
+
+   IMPLICIT NONE
+   PRIVATE
+
+   INTEGER, PUBLIC, PARAMETER ::   jpbyt   = 8    !: real size for mpp communications
+   INTEGER, PUBLIC, PARAMETER ::   jpbytda = 4    !: real size in input data files 4 or 8
+
+   ! Number model from which the SELECTED_*_KIND are requested:
+   !             4 byte REAL       8 byte REAL
+   ! CRAY:           -            precision = 13
+   !                              exponent = 2465
+   ! IEEE:      precision = 6     precision = 15
+   !            exponent = 37     exponent = 307
+
+   !                                                                !!** Floating point **
+   INTEGER, PUBLIC, PARAMETER ::   sp = SELECTED_REAL_KIND( 6, 37)   !: single precision (real 4)
+   INTEGER, PUBLIC, PARAMETER ::   dp = SELECTED_REAL_KIND(12,307)   !: double precision (real 8)
+   INTEGER, PUBLIC, PARAMETER ::   wp = dp                              !: working precision
+
+   !                                                                !!** Integer **
+   INTEGER, PUBLIC, PARAMETER ::   i4 = SELECTED_INT_KIND( 9)        !: single precision (integer 4)
+   INTEGER, PUBLIC, PARAMETER ::   i8 = SELECTED_INT_KIND(14)        !: double precision (integer 8)
+   
+   !                                                                !!** Integer **
+   INTEGER, PUBLIC, PARAMETER ::   lc = 256                          !: Lenght of Character strings
+
+   !!----------------------------------------------------------------------
+   !! NEMO 3.3 , NEMO Consortium (2010)
+   !! $Id$ 
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+END MODULE par_kind
Index: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/TOOLS/RK_PARTITION/src/partition_mod.f90
===================================================================
--- branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/TOOLS/RK_PARTITION/src/partition_mod.f90	(revision 3849)
+++ branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/TOOLS/RK_PARTITION/src/partition_mod.f90	(revision 3849)
@@ -0,0 +1,2264 @@
+MODULE partition_mod
+   USE par_kind
+   IMPLICIT none
+
+   PRIVATE
+
+   LOGICAL, PARAMETER :: PARTIT_DEBUG = .FALSE. ! Ctrl of debug messages from
+                                                ! partitioning code
+   LOGICAL, PARAMETER :: TRIM_DEBUG = .FALSE. ! Ctrl of debug messages from
+                                              ! trimming code.
+
+   INTEGER, PUBLIC, SAVE :: jpni, jpnj, jpnij
+   INTEGER, PUBLIC, SAVE :: jpiglo, jpjglo
+
+  ! Value representing land in the mask used for partitioning and message
+  ! trimming
+  INTEGER, PUBLIC, PARAMETER :: LAND = 0
+
+  INTEGER, PARAMETER :: MAX_FACTORS = 20 ! Max no. of factors factor() can return
+  ! Most of the following parameters should be input vars - TODO ARPDEBUG
+  INTEGER, PARAMETER :: jpreci = 1, jprecj = 1
+  INTEGER, PARAMETER :: nextra = 2 ! Extent to back off from coastlines when
+                                   ! trimming - causes sea-ice problems if 0.
+  LOGICAL, PARAMETER :: trim_land = .TRUE. ! Whether to trim land from sub domains
+  LOGICAL, PUBLIC, SAVE :: cyclic_bc ! Whether we have cyclic bound. cond.'s
+
+  INTEGER, SAVE :: nprocp ! No. of PEs to partition over
+
+   ! Information held by the ensemble leader about all processes.
+   !     pielb      Lower (west) longitude bound index.
+   !     pieub      Upper (east) longitude bound index.
+   !     piesub     Number of longitude gridpoints.
+   !     pilbext    True if the lower longitude boundary is external.
+   !     piubext    True if the upper longitude boundary is external.
+   !     pjelb      Lower (south) latitude bound index.
+   !     pjeub      Upper (north) latitude bound index.
+   !     pjesub     Number of latitude gridpoints.
+   !     pjlbext    True if the lower latitude boundary is external.
+   !     pjubext    True if the upper latitude boundary is external.
+   
+   INTEGER, PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) :: pielb, pieub, piesub, pjelb, pjeub, pjesub, pnactive
+   LOGICAL, PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:) :: pilbext, piubext, pjlbext, pjubext
+
+   ! Process dependent partitioning information.
+
+   !     ielb       Lower (west) longitude bound index.
+   !     ieub       Upper (east) longitude bound index.
+   !     iesub      Number of longitude gridpoints.
+   !     ilbext     True if the lower longitude boundary is external.
+   !     iubext     True if the upper longitude boundary is external.
+   !     jelb       Lower (south) latitude bound index.
+   !     jeub       Upper (north) latitude bound index.
+   !     jesub      Number of latitude gridpoints.
+   !     jlbext     True if the lower latitude boundary is external.
+   !     jubext     True if the upper latitude boundary is external.
+   !     trimmed    True if boundary on given proc trimmed.
+
+   INTEGER, SAVE :: ielb, ieub, iesub
+   INTEGER, SAVE :: jelb, jeub, jesub
+   LOGICAL, SAVE :: ilbext, iubext, jlbext, jubext
+
+   ! Stores whether a domain's boundaries have been trimmed as
+   ! trimmed(boundary, PE) where boundary is one of {n,e,s,w}idx
+   ! for the Northern, Eastern, Southern or Western boundary, respectively.
+   LOGICAL, SAVE, ALLOCATABLE, DIMENSION(:,:) :: trimmed
+   INTEGER, PARAMETER :: nidx = 1, eidx= 2, sidx = 3, widx = 4
+
+   ! Parameters for the cost function used when evaluating different 
+   ! possible domain partitions
+
+   !     Mnemonics:
+   !      p = process (i.e. core)
+   !      n = node
+   !      l = length (number of points)
+   !      c = count  (number of messages)
+   !      i = internal (intra-node, or on-node)
+   !      x = external (inter-node, or off-node)
+
+   INTEGER, PARAMETER :: pv_index_overall = 1
+   INTEGER, PARAMETER :: pv_index_wet     = 2
+   INTEGER, PARAMETER :: pv_index_dry     = 3
+   INTEGER, PARAMETER :: pv_index_pli     = 4
+   INTEGER, PARAMETER :: pv_index_plx     = 5
+   INTEGER, PARAMETER :: pv_index_pci     = 6
+   INTEGER, PARAMETER :: pv_index_pcx     = 7
+   INTEGER, PARAMETER :: pv_index_nli     = 8
+   INTEGER, PARAMETER :: pv_index_nlx     = 9
+   INTEGER, PARAMETER :: pv_index_nci     = 10
+   INTEGER, PARAMETER :: pv_index_ncx     = 11
+   INTEGER, PARAMETER :: pv_index_tli     = 12
+   INTEGER, PARAMETER :: pv_index_tlx     = 13
+   INTEGER, PARAMETER :: pv_index_tci     = 14
+   INTEGER, PARAMETER :: pv_index_tcx     = 15
+
+   INTEGER, PARAMETER :: pv_index_pcomms  = 16
+   INTEGER, PARAMETER :: pv_index_ncomms  = 17
+   INTEGER, PARAMETER :: pv_index_tcomms  = 18
+     
+   INTEGER, PARAMETER :: pv_num_scores    = 18
+   
+   REAL(wp),PARAMETER :: pv_awful = 1.0e20
+
+   PUBLIC allocate_arrays, nemo_partition, write_partition, partition_mca_rk
+
+CONTAINS
+
+   SUBROUTINE allocate_arrays(num_pes, ierr)
+     !!----------------------------------------------------------------------
+     ! Allocate partitioning arrays.
+     !!----------------------------------------------------------------------
+     IMPLICIT none
+     INTEGER, INTENT(In)  :: num_pes
+     INTEGER, INTENT(out) :: ierr
+     !!----------------------------------------------------------------------
+
+      IF ( .NOT. ALLOCATED(pielb) ) THEN
+         ALLOCATE (pielb(num_pes),   pieub(num_pes), piesub(num_pes),     &
+                   pilbext(num_pes), piubext(num_pes),                    &
+                   pjelb(num_pes),   pjeub(num_pes), pjesub(num_pes),     &
+                   pjlbext(num_pes), pjubext(num_pes), pnactive(num_pes), &
+                   trimmed(4,num_pes), Stat = ierr)
+         IF(ierr /= 0)THEN
+            WRITE(*,*)('ERROR: failed to allocate partitioning arrays')
+            RETURN
+         END IF
+      ENDIF
+
+   END SUBROUTINE allocate_arrays
+
+
+   SUBROUTINE nemo_partition( num_pes )
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE nemo_partition  ***
+      !!
+      !! ** Purpose : Work out a sensible factorisation of the number of
+      !!              processors for the x and y dimensions.
+      !! ** Method  :
+      !!----------------------------------------------------------------------
+      IMPLICIT None
+      INTEGER, INTENT(in) :: num_pes ! The number of MPI processes we have
+      !
+      INTEGER             :: ii, jj, jn, ierr
+      INTEGER             :: ifact1, ifact2 ! factors of num_pes, ifact1 <= ifact2
+      INTEGER, ALLOCATABLE, DIMENSION(:,:) :: ilcit, ilcjt, iimppt, ijmppt
+      INTEGER :: nreci, nrecj, iresti, irestj
+      INTEGER :: jpi, jpj, nbondi, nbondj
+      INTEGER :: nldi, nldj, nlei, nlej, nlci, nlcj
+      INTEGER :: zidom, zjdom 
+      !!----------------------------------------------------------------------
+
+      ! Copy no. of PEs into module variable
+      nprocp = num_pes
+
+      ! Factorise the number of processors into ifact1*ifact2, such that
+      ! ifact1 and ifact2 are as nearly equal as possible.
+
+      CALL sqfact( num_pes, ifact1, ifact2 )
+
+      ! Make sure that the smaller dimension of the processor grid
+      ! is given the smaller dimension of the global domain
+      IF( jpiglo <= jpjglo) THEN
+         jpni = ifact1
+         jpnj = ifact2
+      ELSE
+         jpni = ifact2
+         jpnj = ifact1
+      ENDIF
+
+      ! Calc jp{i,j} as done in nemo_init
+      jpi = ( jpiglo-2*jpreci + (jpni-1) ) / jpni + 2*jpreci   ! first  dim.
+      jpj = ( jpjglo-2*jprecj + (jpnj-1) ) / jpnj + 2*jprecj   ! second dim.
+
+      ! This should never happen
+      IF( (jpni*jpnj) /= num_pes) THEN
+         WRITE (*, *) 'ERROR: internal error - factorisation of number of PEs failed'
+         RETURN
+      ENDIF
+
+      ! This should only happen if num_pes is prime
+      IF( ifact1 <= 1 ) THEN
+         WRITE (*, *) 'ERROR: factorisation of number of PEs failed'
+         RETURN
+      ENDIF
+      !
+      jpnij = jpni*jpnj
+
+      WRITE(*,FMT="(/'Generating REGULAR partition on ',I4,'x',I4,' process grid'/)") &
+             jpni, jpnj
+      !
+      !  1. Dimension arrays for subdomains
+      ! -----------------------------------
+      !  Computation of local domain sizes ilcit() ilcjt()
+      !  These dimensions depend on global sizes jpni,jpnj and jpiglo,jpjglo
+      !  The subdomains are squares lesser than or equal to the global
+      !  dimensions divided by the number of processors minus the overlap
+      !  array (cf. par_oce.F90).
+      ALLOCATE(ilcit(jpni,jpnj), ilcjt(jpni,jpnj), &
+               iimppt(jpni,jpnj), ijmppt(jpni,jpnj), Stat=ierr)
+
+      nreci  = 2 * jpreci
+      nrecj  = 2 * jprecj
+      iresti = MOD( jpiglo - nreci , jpni )
+      irestj = MOD( jpjglo - nrecj , jpnj )
+
+      IF(  iresti == 0 )   iresti = jpni
+      DO jj = 1, jpnj
+         DO ii = 1, iresti
+            ilcit(ii,jj) = jpi
+         END DO
+         DO ii = iresti+1, jpni
+            ilcit(ii,jj) = jpi -1
+         END DO
+      END DO
+      
+      IF( irestj == 0 )   irestj = jpnj
+      DO ii = 1, jpni
+         DO jj = 1, irestj
+            ilcjt(ii,jj) = jpj
+         END DO
+         DO jj = irestj+1, jpnj
+            ilcjt(ii,jj) = jpj -1
+         END DO
+      END DO
+ 
+      zidom = nreci
+      DO ii = 1, jpni
+         zidom = zidom + ilcit(ii,1) - nreci
+      END DO
+      !WRITE(*,*)
+      !WRITE(*,*)' sum ilcit(i,1) = ', zidom, ' jpiglo = ', jpiglo
+      
+      zjdom = nrecj
+      DO jj = 1, jpnj
+         zjdom = zjdom + ilcjt(1,jj) - nrecj
+      END DO
+      !WRITE(*,*)' sum ilcit(1,j) = ', zjdom, ' jpjglo = ', jpjglo
+      !WRITE(*,*)
+      
+
+      !  2. Index arrays for subdomains
+      ! -------------------------------
+      
+      iimppt(:,:) = 1
+      ijmppt(:,:) = 1
+      
+      IF( jpni > 1 ) THEN
+         DO jj = 1, jpnj
+            DO ii = 2, jpni
+               iimppt(ii,jj) = iimppt(ii-1,jj) + ilcit(ii-1,jj) - nreci
+            END DO
+         END DO
+      ENDIF
+
+      IF( jpnj > 1 ) THEN
+         DO jj = 2, jpnj
+            DO ii = 1, jpni
+               ijmppt(ii,jj) = ijmppt(ii,jj-1)+ilcjt(ii,jj-1)-nrecj
+            END DO
+         END DO
+      ENDIF
+      
+      ! 3. Subdomain description
+      ! ------------------------
+
+      DO jn = 1, jpnij
+         ii = 1 + MOD( jn-1, jpni )
+         jj = 1 + (jn-1) / jpni
+         nlci       =  ilcit (ii,jj)  
+         nlcj       =  ilcjt (ii,jj)   
+
+         pielb(jn) = iimppt(ii,jj)
+         pjelb(jn) = ijmppt(ii,jj)
+
+         nbondj = -1                                   ! general case
+         IF( jn   >  jpni          )   nbondj = 0      ! first row of processor
+         IF( jn   >  (jpnj-1)*jpni )   nbondj = 1      ! last  row of processor
+         IF( jpnj == 1             )   nbondj = 2      ! one processor only in j-direction
+         nbondi = 0                                    ! 
+         IF( MOD( jn, jpni ) == 1 )   nbondi = -1      !
+         IF( MOD( jn, jpni ) == 0 )   nbondi =  1      !
+         IF( jpni            == 1 )   nbondi =  2      ! one processor only in i-direction
+         nldi =  1   + jpreci
+         nlei = nlci - jpreci
+         pilbext(jn) = .FALSE.
+         piubext(jn) = .FALSE.
+         IF( nbondi == -1 .OR. nbondi == 2 )THEN
+            nldi = 1
+            pilbext(jn) = .TRUE.
+         END IF
+         IF( nbondi ==  1 .OR. nbondi == 2 )THEN
+            nlei = nlci
+            piubext(jn) = .TRUE.
+         END IF
+
+         nldj =  1   + jprecj
+         nlej = nlcj - jprecj
+         pjlbext(jn) = .FALSE.
+         pjubext(jn) = .FALSE.
+         IF( nbondj == -1 .OR. nbondj == 2 )THEN
+            nldj = 1
+            pjlbext(jn) = .TRUE.
+         END IF
+         IF( nbondj ==  1 .OR. nbondj == 2 )THEN
+            nlej = nlcj
+            pjubext(jn) = .TRUE.
+         END IF
+
+         ! nldj corresponds to domain lower bound and nlej the upper bound
+         pieub(jn) = pielb(jn) + (nlei - nldi + 1) - 1
+         pjeub(jn) = pjelb(jn) + (nlej - nldj + 1) - 1
+
+         ! No trimming in this partitioning
+         trimmed(:,jn) = .FALSE.
+      END DO
+
+   END SUBROUTINE nemo_partition
+
+   !==========================================================================
+
+   SUBROUTINE sqfact ( kn, kna, knb )
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE sqfact  ***
+      !!
+      !! ** Purpose :   return factors (kna, knb) of kn, such that
+      !!                (1) kna*knb=kn
+      !!                (2) kna and knb are as near equal as possible
+      !!                (3) kna < knb
+      !! ** Method  :   Search backwards from the square root of kn,
+      !!                until we find an integer that cleanly divides kn
+      !! ** Preconditions : kn must be positive
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in   ) ::   kn
+      INTEGER, INTENT(  out) ::   kna, knb
+       
+      ! Search backwards from the square root of n.
+
+      fact_loop: DO kna=INT(SQRT(REAL(kn))),1,-1
+         IF ( kn/kna*kna == kn ) THEN
+            EXIT fact_loop
+         ENDIF
+      END DO fact_loop
+
+      IF( kna < 1 ) kna = 1 
+
+      ! kna divides kn cleanly. Work out the other factor.
+      knb = kn/kna
+
+   END SUBROUTINE sqfact
+
+   !==========================================================================
+
+   SUBROUTINE partition_mca_rk(mask, nproc, ncpn, ierr)
+     IMPLICIT NONE
+     !!------------------------------------------------------------------
+     !! Multi-Core Aware recursive partitioning of the domain. As for 
+     !! partition_rk but choose the partion 
+     !! so as to minimize off-node MPI communication
+     !!
+     !! Original by Stephen Pickles for POLCOMS code.
+     !! Implementation in NEMO by Andrew Porter, 26/01/2012
+     !!------------------------------------------------------------------
+     ! Subroutine arguments.
+     INTEGER, INTENT(in)        :: mask(:,:)
+     INTEGER, INTENT(in)        :: nproc
+     INTEGER, INTENT(in)        :: ncpn ! No. of cores per node on target machine
+     INTEGER, INTENT(out)       :: ierr
+     ! Local variables
+     INTEGER :: ii
+     INTEGER, DIMENSION(MAX_FACTORS) :: fx, fy, factors
+     INTEGER, DIMENSION(MAX_FACTORS) :: df, multiplicity
+     INTEGER :: nfx, nfy, nfactors, ndf, nperms
+     !INTEGER :: check_nprocx, check_nprocy, check_nprocp
+     INTEGER :: iperm
+     CHARACTER(LEN=256) :: fstr
+     INTEGER :: nprocx, nprocy
+
+     ! A high score is bad
+     REAL(wp), DIMENSION(pv_num_scores)   :: score
+     REAL(wp) :: best_score
+     INTEGER  :: best_perm
+     REAL(wp), DIMENSION(2,pv_num_scores) :: best, wrst
+     !!------------------------------------------------------------------
+
+     ! Copy no. of PEs into module variable
+     nprocp = nproc
+
+     ! Factorise the total number of MPI processes that we have
+     CALL factor (factors,MAX_FACTORS,nfactors,nproc,ierr)
+
+     IF ( ierr.NE.0 ) THEN
+        WRITE (*,FMT="(/'partition_mca_rk: factorisation failed ',I3)") ierr
+        RETURN
+     ELSE
+        WRITE (*,FMT="(/'partition_mca_rk: factors are: ',(I3,1x))") factors(:nfactors)
+        WRITE (*,*)
+     ENDIF
+
+     CALL calc_perms( nfactors, factors,     &
+                      ndf, df, multiplicity, &
+                      nperms )
+
+     DO ii=1,pv_num_scores
+        best(1,ii) = pv_awful
+        best(2,ii) = -1.0_wp
+     END DO
+     DO ii=1,pv_num_scores
+        wrst(1,ii) = 0.0_wp
+        wrst(2,ii) = -1.0_wp
+     END DO
+
+     WRITE(*,"('% Partn',2X,10(A4,2X),4(A9,1X),A7)")       &
+          'Wet', 'Dry',                               &
+          'pli', 'plx', 'pci', 'pcx',                 &
+          'nlx', 'ncx', 'tlx', 'tcx',                 &
+          'P comms', 'N comms', 'T comms', 'Overall', &
+          'Factors'
+
+     !       perm_loop: DO iperm=myinst, nperms-1, nprocp
+     ! Now serial so we have to do them all ourselves - could use
+     ! OpenMP?
+     perm_loop: DO iperm=1, nperms-1, 1
+
+        CALL get_perm_factors( iperm, nfactors, ndf, df, multiplicity, &
+                               fx, nfx, fy, nfy,                       &
+                               nprocx, nprocy)
+
+        CALL partition_rk_core(mask, jpiglo, jpjglo,    &
+             MAX_FACTORS, fx, nfx, fy, nfy, nprocp, ierr)
+
+        IF (ierr .NE. 0) CYCLE perm_loop
+        CALL finish_partition(mask)
+
+        ! Compute the cost function for this partition
+        !
+        CALL eval_partition(ncpn, jpiglo, jpjglo, mask, score)
+        CALL factor_string(fstr,nfx,fx,nfy,fy)
+
+        WRITE (6,'(''%'',I6,1X,10(I5,1X),3(F9.2,1X),E12.4,1x,(A))') &
+             iperm,                                     &
+             INT(score(pv_index_wet)),                  &
+             INT(score(pv_index_dry)),                  &
+             INT(score(pv_index_pli)),                  &
+             INT(score(pv_index_plx)),                  &
+             INT(score(pv_index_pci)),                  &
+             INT(score(pv_index_pcx)),                  &
+             INT(score(pv_index_nlx)),                  &
+             INT(score(pv_index_ncx)),                  &
+             INT(score(pv_index_tlx)),                  &
+             INT(score(pv_index_tcx)),                  &
+             score(pv_index_pcomms),                    &
+             score(pv_index_ncomms),                    &
+             score(pv_index_tcomms),                    &
+             score(pv_index_overall),                   &
+             TRIM(fstr)
+
+        DO ii=1,pv_num_scores
+           IF (score(ii) .LT.  best(1,ii)) THEN
+              best(1,ii) = score(ii)
+              best(2,ii) = iperm
+           END IF
+           IF (score(ii) .GT. wrst(1,ii)) THEN
+              wrst(1,ii) = score(ii)
+              wrst(2,ii) = iperm
+           END IF
+        END DO
+
+     END DO perm_loop
+
+     !  Now choose the "best" partition
+
+     best_score = best(1,pv_index_overall)
+     best_perm  = INT( best(2,pv_index_overall) )
+
+     WRITE (*,'(/A32,A20,A20)')  &
+          ' ','  best score / perm ','  worst score / perm'
+     WRITE (*,'(A32,F12.3,F8.0,E12.3,F8.0)') 'overall: ',    &
+          best(1,pv_index_overall), best(2,pv_index_overall), &
+          wrst(1,pv_index_overall), wrst(2,pv_index_overall)
+     WRITE (*,'(A32,F12.3,F8.0,E12.3,F8.0)') 'wet points: ', &
+          best(1,pv_index_wet), best(2,pv_index_wet),      &
+          wrst(1,pv_index_wet), wrst(2,pv_index_wet)
+     WRITE (*,'(A32,F12.3,F8.0,E12.3,F8.0)') 'dry points: ', &
+          best(1,pv_index_dry), best(2,pv_index_dry),      &
+          wrst(1,pv_index_dry), wrst(2,pv_index_dry)
+     WRITE (*,'(A32,F12.3,F8.0,E12.3,F8.0)')                 &
+          'proc max  on-node wet points: ',                &
+          best(1,pv_index_pli), best(2,pv_index_pli),      &
+          wrst(1,pv_index_pli), wrst(2,pv_index_pli)
+     WRITE (*,'(A32,F12.3,F8.0,E12.3,F8.0)')                 &
+          'proc max off-node wet points: ',                &
+          best(1,pv_index_plx), best(2,pv_index_plx),      &
+          wrst(1,pv_index_plx), wrst(2,pv_index_plx)
+     WRITE (*,'(A32,F12.3,F8.0,E12.3,F8.0)')                 &
+          'proc max  on-node   messages: ',                &
+          best(1,pv_index_pci), best(2,pv_index_pci),      &
+          wrst(1,pv_index_pci), wrst(2,pv_index_pci)
+     WRITE (*,'(A32,F12.3,F8.0,E12.3,F8.0)')                 &
+          'proc max off-node   messages: ',                &
+          best(1,pv_index_pcx), best(2,pv_index_pcx),      &
+          wrst(1,pv_index_pcx), wrst(2,pv_index_pcx)
+     WRITE (*,'(A32,F12.3,F8.0,E12.3,F8.0)')                 &
+          'node max off-node wet points: ',                &
+          best(1,pv_index_nlx), best(2,pv_index_nlx),      &
+          wrst(1,pv_index_nlx), wrst(2,pv_index_nlx)
+     WRITE (*,'(A32,F12.3,F8.0,E12.3,F8.0)')                 &
+          'node max off-node   messages: ',                     &
+          best(1,pv_index_ncx), best(2,pv_index_ncx),         &
+          wrst(1,pv_index_ncx), wrst(2,pv_index_ncx)
+     WRITE (*,'(A32,F12.3,F8.0,E12.3,F8.0)')                 &
+          'total    off-node wet points: ',                     &
+          best(1,pv_index_tlx), best(2,pv_index_tlx),         &
+          wrst(1,pv_index_tlx), wrst(2,pv_index_tlx)
+     WRITE (*,'(A32,F12.3,F8.0,E12.3,F8.0)')                 &
+          'per core communications cost: ',                     &
+          best(1,pv_index_pcomms), best(2,pv_index_pcomms),   &
+          wrst(1,pv_index_pcomms), wrst(2,pv_index_pcomms)
+     WRITE (*,'(A32,F12.3,F8.0,E12.3,F8.0)')                 &
+          'per node communications cost: ',                     &
+          best(1,pv_index_ncomms), best(2,pv_index_ncomms),   &
+          wrst(1,pv_index_ncomms), wrst(2,pv_index_ncomms)
+     WRITE (*,'(A32,F12.3,F8.0,E12.3,F8.0)')                 &
+          'network communications cost: ',                      &
+          best(1,pv_index_tcomms), best(2,pv_index_tcomms),   &
+          wrst(1,pv_index_tcomms), wrst(2,pv_index_tcomms)
+
+     WRITE (*,"(/'partition_mca_rk: overall best perm is ',I6,', score=',F12.3/)") &
+          best_perm, best_score
+
+     ! Use the same partition on all processes
+
+     ! If a particular factorisation has been forced by
+     ! the nn_{x,y}factors fields in the nammpp section of the namelist
+     ! then use that one instead
+
+!!$      IF ((nxfactors + nyfactors) > 0) THEN
+!!$         check_nprocx = 1
+!!$         check_nprocy = 1
+!!$         DO ii=1,nxfactors
+!!$            check_nprocx = check_nprocx*xfactors(ii)
+!!$         END DO
+!!$         DO ii=1,nyfactors
+!!$            check_nprocy = check_nprocy*yfactors(ii)
+!!$         END DO
+!!$         check_nprocp = check_nprocx*check_nprocy
+!!$         IF (check_nprocp .EQ. nprocp) THEN
+!!$            nprocx = check_nprocx
+!!$            nprocy = check_nprocy
+!!$            nfx = nxfactors
+!!$            nfy = nyfactors
+!!$            fx(1:nfx) = xfactors(1:nfx)
+!!$            fy(1:nfy) = yfactors(1:nfy)
+!!$         ELSE
+!!$            CALL ctl_stop('STOP', 'part_mca_rk: '//   &
+!!$                          'requested factorisation does not match nprocp')
+!!$         END IF
+!!$      ELSE
+     ! Use the best factorisation found above
+     IF (best_perm < 0.0_wp) THEN
+        WRITE (*,*) 'partition_mca_rk: no feasible partition found'
+        ierr = 1
+        RETURN
+     END IF
+     CALL get_perm_factors(best_perm, nfactors, ndf, df, multiplicity, &
+          fx, nfx, fy, nfy,                           &
+          nprocx, nprocy)
+!!$      END IF
+
+     ! Set corresponding NEMO variables for PE grid, even though it is now
+     ! rather irregular
+     jpni = nprocx
+     jpnj = nprocy
+
+     WRITE (*,*)
+     WRITE (*,'(A39)',advance='no') &
+          'partition_mca_rk: partitioning with factors '
+     CALL print_factors(nfx, fx, nfy, fy)
+
+     CALL partition_rk_core ( mask, jpiglo, jpjglo, &
+          MAX_FACTORS,          &
+          fx, nfx, fy, nfy,     &
+          nproc, ierr )
+
+     IF (ierr .NE. 0) THEN
+        WRITE (*,*) 'partition_mca_rk: failed to apply selection partition'
+        RETURN
+     END IF
+     CALL finish_partition(mask)
+
+     CALL eval_partition(ncpn, jpiglo, jpjglo, mask, score)
+     CALL factor_string(fstr,nfx,fx,nfy,fy)
+     WRITE(*,'(10(A7,1X),4(A9,1X),A7)')           &
+          'Wet', 'Dry',                               &
+          'pli', 'plx', 'pci', 'pcx',                 &
+          'nlx', 'ncx', 'tlx', 'tcx',                 &
+          'P comms', 'N comms', 'T comms', 'Overall', &
+          'Factors'
+     WRITE (*,'(10(F7.0,1X),4(F9.2,1X),A50)') &
+          score(pv_index_wet),              &
+          score(pv_index_dry),              &
+          score(pv_index_pli),              &
+          score(pv_index_plx),              &
+          score(pv_index_pci),              &
+          score(pv_index_pcx),              &
+          score(pv_index_nlx),              &
+          score(pv_index_ncx),              &
+          score(pv_index_tlx),              &
+          score(pv_index_tcx),              &
+          score(pv_index_pcomms),           &
+          score(pv_index_ncomms),           &
+          score(pv_index_tcomms),           &
+          score(pv_index_overall),          &
+          fstr
+
+   END SUBROUTINE partition_mca_rk
+
+
+   SUBROUTINE partition_rk_core( mask, nx, ny, maxfax, fx, nfx, fy, nfy,   &
+        nproc, ierr )
+     IMPLICIT NONE
+     !!------------------------------------------------------------------
+     !!------------------------------------------------------------------
+     INTEGER, INTENT(in)  :: nx, ny
+     INTEGER, INTENT(in)  :: mask(nx,ny)
+     INTEGER, INTENT(in)  :: maxfax, nfx, nfy
+     INTEGER, INTENT(in)  :: fx(maxfax), fy(maxfax)
+     INTEGER, INTENT(in)  :: nproc
+     INTEGER, INTENT(out) :: ierr
+     ! Local variables
+     INTEGER :: istart, jstart, istop, jstop
+     INTEGER :: f,gnactive 
+     INTEGER :: i,ifax,nfax,ifin,ifx,ify,ilb,iproc
+     INTEGER :: ist,isub,isub_old,isub_new,iub    
+     INTEGER :: j,jfin,jlb,jst,jub,line
+     INTEGER :: ngone,nsub_old,nsub_new,ntarget,ntry
+     LOGICAL :: partx
+
+     ! Zero the error flag
+     ierr = 0
+
+     ! Count the significant (non-zero) factors.
+     nfax = nfx+nfy
+     IF( PARTIT_DEBUG )THEN
+        WRITE (*,"('jpni = ',I3,1x,'nfx = ',I3,/'Factorized into: ',(I3,1x))") &
+             jpni, nfx, fx(:nfx)
+        WRITE (*,"('jpnj = ',I3,1x,'nfy = ',I3,/'Factorized into: ',(I3,1x))") & 
+             jpnj, nfy, fy(:nfy)
+        WRITE (*,"('Partitioning a total of ',I3,' times')") nfax
+     END IF
+
+     ! Define the full domain as the one and only sub-domain.
+     istart = 1
+     istop = nx
+     jstart = 1
+     jstop = ny
+
+     nsub_old = 1
+     pielb(nsub_old) = istart
+     pjelb(nsub_old) = jstart
+     pieub(nsub_old) = istop
+     pjeub(nsub_old) = jstop
+     ! Count the number of active points.
+     gnactive = 0
+     DO j=jstart,jstop,1
+        DO i=istart,istop,1
+           IF ( mask(i,j) == 1 ) gnactive = gnactive+1
+        ENDDO
+     ENDDO
+
+     IF( PARTIT_DEBUG )THEN
+       WRITE (*,*) 'Total wet points ',gnactive
+     END IF
+
+     pnactive(nsub_old) = gnactive
+
+     ! Partition for each factor.
+     DO ifax=1,nfax
+        IF ( ifax.EQ.1 ) THEN
+           ! Start by partitioning the dimension with more factors.
+           partx = nfx.GE.nfy
+           ifx = 0
+           ify = 0
+        ELSE
+           ! Try to toggle the partitioning direction.
+           partx = .NOT.partx
+           IF ( partx ) THEN
+              ! If we have run out of factors in x, switch to y.
+              partx = .NOT. ifx+1.GT.nfx
+           ELSE
+              ! If we have run out of factors in y, switch to x.
+              partx =       ify+1.GT.nfy
+           ENDIF
+        ENDIF
+
+        IF( PARTIT_DEBUG )THEN
+          WRITE (*,*)
+          WRITE (*,*) '###########################################'
+          WRITE (*,*)
+          WRITE (*,*) 'Partition ',ifax,'partx = ',partx
+        END IF
+
+        IF ( partx ) THEN
+           ! ============================================================
+           !         Partition in x.
+           ! ============================================================
+           ifx = ifx+1
+           f = fx(ifx)
+           nsub_new = nsub_old*f
+
+           IF( PARTIT_DEBUG )THEN
+              WRITE (*,*) 'Partitioning in x from ',nsub_old &
+                   ,' to ',nsub_new
+           END IF
+
+           DO isub_old=1, nproc, nproc/nsub_old
+              ! Check that there are sufficient points to factorise.
+              IF ( pieub(isub_old)-pielb(isub_old)+1.LT.f ) THEN
+                 WRITE (*,*)  &
+                      'partition_rk: insufficient points to partition'
+                 ierr = 1
+                 RETURN
+              ENDIF
+
+              ! Set the target number of points in the new sub-domains.
+              ntarget = NINT(REAL(pnactive(isub_old))/REAL(f))
+              ist = pielb(isub_old)
+              iub = pieub(isub_old)
+              jlb = pjelb(isub_old)
+              jub = pjeub(isub_old)
+
+              IF( PARTIT_DEBUG )THEN
+                WRITE (*,"(/'=======================================')")
+                WRITE (*,"(/'Partitioning sub-domain ',I3,': (',I3,':',I3,')(',I3,':',I3,')')") &
+                     isub_old,ist,iub,jlb,jub
+                WRITE (*,"('Old partition has ',I8,' points')") pnactive(isub_old)
+                WRITE (*,"('Target is ',I8,' points')") ntarget
+              END IF
+              ! Create the new sub-domains.
+              ngone = 0
+              DO isub=1,f,1
+                 isub_new = isub_old + (isub-1)*nproc/nsub_new
+
+                 IF( PARTIT_DEBUG )THEN
+                    WRITE (*,*)
+                    WRITE (*,*) 'Creating new sub-domain ',isub_new
+                 END IF
+                 ! The last new sub-domain takes the remaining data.
+                 IF ( isub.EQ.f ) THEN
+                    ifin = iub
+                 ELSE
+                    ! Scan lines in x counting active grid points until we
+                    ! exceed the target.
+                    ntry = 0
+                    ifin = -1
+                    scanrows: DO i=ist,iub,1
+                       ! Count up the contribution from the next line.
+                       line = 0
+                       DO j=jlb,jub,1
+                          IF ( mask(i,j) == 1 ) line = line+1
+                       ENDDO
+
+                       ! Check against the target.
+                       IF ( ntry+line.GT.ntarget ) THEN
+                          ! Is it nearer the target to include this line or not ?
+                          IF ( ntry+line-ntarget.GT.ntarget-ntry ) THEN
+                             ! Nearer start than end. Finish at previous line.
+                             ifin = i-1
+                          ELSE
+                             ! Nearer end than start. Include this line.
+                             ifin = i
+                             ntry = ntry + line
+                          ENDIF
+
+                          IF( PARTIT_DEBUG )THEN
+                             WRITE (*,*) 'Reached target at ',ifin  &
+                                  ,' ntry = ',ntry
+                          END IF
+                          EXIT scanrows
+                       ENDIF
+                       ! Add in the line count to the running total.
+                       ntry = ntry + line
+                    ENDDO scanrows
+                    IF ( ifin.LT.0 ) ifin = iub
+                 ENDIF
+
+                 ! Set up the parameters of the new sub-domain.
+                 pnactive(isub_new) = 0
+                 DO j=jlb,jub
+                    DO i=ist,ifin
+                       IF ( mask(i,j) == 1 ) THEN
+                          pnactive(isub_new) = pnactive(isub_new)+1
+                       ENDIF
+                    ENDDO
+                 ENDDO
+                 pielb(isub_new) = ist
+                 pieub(isub_new) = ifin
+                 pjelb(isub_new) = jlb
+                 pjeub(isub_new) = jub
+
+                 IF( PARTIT_DEBUG )THEN
+                    WRITE (*,*) 'New subdomain is ',ist,ifin,jlb,jub &
+                         ,' with ',pnactive(isub_new),' points'
+                 END IF
+                 !             Reset the target based on the actual number of points
+                 !             remaining, split between the partitions remaining.
+                 ngone = ngone+ntry
+                 !             Start the next search at the next point.
+                 ist = ifin+1
+              ENDDO
+           ENDDO
+
+        ELSE
+
+           !      ============================================================
+           !      Partition in y.
+           !      ============================================================
+           ify = ify+1
+           f = fy(ify)
+           nsub_new = nsub_old*f
+
+           IF( PARTIT_DEBUG )THEN
+              WRITE (*,*) 'Partitioning in y from ',nsub_old &
+                   ,' to ',nsub_new
+           END IF
+
+           DO isub_old=1, nproc, nproc/nsub_old
+              !           Check that there are sufficient points to factorise.
+              IF ( pjeub(isub_old)-pjelb(isub_old)+1.LT.f ) THEN
+                 WRITE (*,*)  &
+                      'partition_rk: insufficient points to partition'
+                 ierr = 1
+                 RETURN
+              ENDIF
+              !           Set the target number of points in the new sub-domains.
+              ntarget = NINT(REAL(pnactive(isub_old))/REAL(f))
+              ilb = pielb(isub_old)
+              iub = pieub(isub_old)
+              jst = pjelb(isub_old)
+              jub = pjeub(isub_old)
+
+              IF( PARTIT_DEBUG )THEN
+                 WRITE (*,*)
+                 WRITE (*,*) '======================================='
+                 WRITE (*,*)
+                 WRITE (*,*) 'Partitioning sub-domain ',isub_old,' : ' &
+                      ,ilb,iub,jst,jub
+                 WRITE (*,*) 'Old partition has ',pnactive(isub_old)   &
+                      ,' points'
+                 WRITE (*,*) 'Target is ',ntarget
+              END IF
+              !           Create the new sub-domains.
+              ngone = 0
+              DO isub=1,f
+                 isub_new = isub_old + (isub-1)*nproc/nsub_new
+
+                 IF( PARTIT_DEBUG )THEN
+                    WRITE (*,*)
+                    WRITE (*,*) 'Creating new sub-domain ',isub_new
+                 END IF
+
+                 !             The last new sub-domain takes the remaining data.
+                 IF ( isub == f ) THEN
+                    jfin = jub
+                 ELSE
+                    !      Scan lines in y counting active grid points until we
+                    !      exceed the target.
+                    ntry = 0
+                    jfin = -1
+                    scancols: DO j=jst,jub
+                       !     Count up the contribution from the next line.
+                       line = 0
+                       DO i=ilb,iub
+                          IF ( mask(i,j) == 1 ) line = line+1
+                       ENDDO
+
+                       ! Check against the target.
+                       IF ( ntry+line > ntarget ) THEN
+                          ! Is it nearer the target to include this line or not ?
+                          IF ( ntry+line-ntarget > ntarget-ntry ) THEN
+                             ! Nearer start than end. Finish at previous line.
+                             jfin = j-1
+                          ELSE
+                             ! Nearer end than start. Include this line.
+                             jfin = j
+                             ntry = ntry + line
+                          ENDIF
+
+                          IF( PARTIT_DEBUG )THEN
+                             WRITE (*,*) 'Reached target at ',jfin &
+                                  ,' ntry = ',ntry
+                          END IF
+
+                          EXIT scancols
+                       ENDIF
+                       ! Add in the line count to the running total.
+                       ntry = ntry + line
+                    ENDDO scancols
+                    IF ( jfin.LT.0 ) jfin = jub
+                 ENDIF
+                 ! Set up the parameters of the new sub-domain.
+                 pnactive(isub_new) = 0
+                 DO j=jst,jfin
+                    DO i=ilb,iub
+                       IF ( mask(i,j) == 1 ) THEN
+                          pnactive(isub_new) = pnactive(isub_new)+1
+                       ENDIF
+                    ENDDO
+                 ENDDO
+                 pielb(isub_new) = ilb
+                 pieub(isub_new) = iub
+                 pjelb(isub_new) = jst
+                 pjeub(isub_new) = jfin
+
+                 IF( PARTIT_DEBUG )THEN
+                    WRITE (*,*) 'New subdomain is ',ilb,iub,jst,jfin &
+                         ,' with ',pnactive(isub_new),' points'
+                 END IF
+
+                 !        Reset the target based on the actual number of points
+                 !        remaining, split between the partitions remaining.
+                 ngone = ngone+ntry
+
+                 IF( PARTIT_DEBUG ) WRITE (*,*) 'Target reset to ',ntarget
+
+                 !             Start the next search at the next point.
+                 jst = jfin+1
+              ENDDO
+           ENDDO
+        ENDIF
+        !       Update the number of sub-domains ready for the next loop.
+        nsub_old = nsub_new
+     ENDDO
+
+     IF( PARTIT_DEBUG )THEN
+        WRITE (*,*)
+        WRITE (*,*) '========================================='
+        WRITE (*,*)
+     END IF
+
+     !     Set the size of each sub-domain.
+     DO iproc=1, nproc, 1
+        piesub(iproc) = pieub(iproc)-pielb(iproc)+1
+        pjesub(iproc) = pjeub(iproc)-pjelb(iproc)+1
+
+        IF( PARTIT_DEBUG )THEN
+           WRITE (*,*) 'Domain ',iproc,' has ',pnactive(iproc), ' ocean points'
+           WRITE(*,"(I3,': [',I3,':',I3,'][',I3,':',I3,']')") &
+                iproc, pielb(iproc), pieub(iproc), pjelb(iproc), pjeub(iproc)
+        END IF
+     ENDDO
+
+   END SUBROUTINE partition_rk_core
+
+
+   SUBROUTINE factor ( ifax, maxfax, nfax, n, ierr )
+
+     !!------------------------------------------------------------------
+     ! Subroutine to return the prime factors of n.
+     ! nfax factors are returned in array ifax which is of maximum
+     ! dimension maxfax.
+     !!------------------------------------------------------------------
+     IMPLICIT NONE
+     ! Subroutine arguments
+     INTEGER ierr, n, nfax, maxfax
+     INTEGER ifax(maxfax)
+     ! Local variables.
+     INTEGER i, ifac, l, nu
+     INTEGER lfax(20)
+     ! lfax contains the set of allowed factors.
+     DATA (lfax(i),i=1,20) / 71,67,59,53,47,43,41,37,31,29 &
+                            ,23,19,17,13,11, 7, 5, 3, 2, 1 /
+     ! Clear the error flag.
+     ierr = 0
+     ! Find the factors of n.
+     IF ( n.EQ.1 ) THEN
+        nfax = 0
+        GOTO 20
+     ENDIF
+     ! nu holds the unfactorised part of the number.
+     ! nfax holds the number of factors found.
+     ! l points to the allowed factor list.
+     ! ifac holds the current factor.
+     nu = n
+     nfax = 0
+     l = 1
+     ifac = lfax(l)
+     ! Label 10 is the start of the factor search loop.
+10   CONTINUE
+     ! Test whether the factor will divide.
+     IF ( MOD(nu,ifac).EQ.0 ) THEN
+        ! Add the factor to the list.
+        nfax = nfax+1
+        IF ( nfax.GT.maxfax ) THEN
+           ierr = 6
+           WRITE (*,*)  &
+                'FACTOR: insufficient space in factor array ',nfax
+           RETURN
+        END IF
+        ifax(nfax) = ifac
+        ! Divide out the factor from the remaining number.
+        ! If unity remains, the number has been
+        ! successfully factorized.
+        nu = nu/ifac
+        IF ( nu.EQ.1   ) GOTO 20
+        ! Loop to try the factor again.
+        GOTO 10
+     ENDIF
+     ! Move on to the next factor in the list.
+     l = l+1
+     ifac = lfax(l)
+     ! Unless the end of the factor list has been reached, loop.
+     IF ( ifac.GT.1 ) GOTO 10
+     ! There is a factor in n which is not in the lfac list.
+     ! Add the remaining number onto the end of the list.
+     nfax = nfax+1
+     IF ( nfax.GT.maxfax ) THEN
+        ierr = 6
+        WRITE (*,*) 'FACTOR: insufficient space in factor array ',nfax
+        RETURN
+     ENDIF
+     ifax(nfax) = nu
+     ! Label 20 is the exit point from the factor search loop.
+20   CONTINUE
+
+   END SUBROUTINE factor
+
+
+   SUBROUTINE part_trim ( depth, isTrimmed, ierr )
+        !!------------------------------------------------------------------
+        !!
+        !! Examines all the sub-domains and trims off boundary rows and
+        !! columns which are all land.
+        !!
+        !!     depth                   real  input     Depth map.
+        !!     isTrimmed               logical output  Whether N,E,S,W edge
+        !!                                             of domain has been 
+        !!                                             trimmed
+        !!     ierr                    int   output    Error flag.
+        !!
+        !!        Mike Ashworth, CLRC Daresbury Laboratory, July 1999
+        !!------------------------------------------------------------------
+        IMPLICIT NONE
+
+        ! Subroutine arguments.
+        INTEGER, INTENT(in)  :: depth(jpiglo,jpjglo)
+        LOGICAL, DIMENSION(:,:), INTENT(out) :: isTrimmed
+        INTEGER, INTENT(out) :: ierr
+        ! Local variables.       
+        INTEGER :: i, iproc, j, newbound
+
+        ! Clear the error code.
+        ierr = 0
+
+        ! Clear all flags
+        ! N E S W
+        ! 1 2 3 4
+        isTrimmed(1:4,1:nprocp) = .FALSE.
+
+        ! Examine each sub-domain in turn.
+        
+        subdomain_loop: DO iproc=1,nprocp
+
+           ! Do not trim if there are no active points at all.
+           ! Otherwise we will trim away the whole sub-domain and we
+           ! will be in big trouble.
+
+           ! Look at the low-i columns (Western edge of domain)
+
+           newbound = pielb(iproc)
+           lowi: DO i=pielb(iproc),pieub(iproc)
+
+              ! Scan this column for wet points
+
+              DO j=MAX(1,pjelb(iproc)-jprecj),MIN(jpjglo,pjeub(iproc)+jprecj)
+
+                 IF ( depth(i,j) == 1 ) THEN
+                     newbound = MAX(i - jpreci - nextra, pielb(iproc))
+
+                     IF( PARTIT_DEBUG )THEN
+                        WRITE (*,*) 'Sub-domain ',iproc,': Low-i loop: row ',i &
+                                     ,' is land. New bound is ',newbound
+                     END IF
+
+                     EXIT lowi
+                 ENDIF
+               ENDDO
+            ENDDO lowi
+
+            ! Update with the new boundary and monitor.
+
+            IF ( newbound.NE.pielb(iproc) ) THEN
+
+               IF( TRIM_DEBUG )THEN
+                  IF ( newbound-pielb(iproc).GT.1 ) THEN
+                     WRITE(*,'(a,i5,3(a,i6))') ' Process ',iproc-1 &
+                                    ,' trimmed ',newbound-pielb(iproc) &
+                          ,' cols: ',pielb(iproc),' to ',newbound-1
+                  ELSE
+                     WRITE(*,'(a,i5,3(a,i6))') ' Process ',iproc-1 &
+                                    ,' trimmed ',newbound-pielb(iproc) &
+                                    ,' col : ',pielb(iproc)
+                  ENDIF
+               END IF
+
+               pielb(iproc) = newbound
+               isTrimmed(4,iproc) = .TRUE.
+            ENDIF
+
+            ! Look at the high-i columns (Eastern edge of domain).
+
+            newbound = pieub(iproc)
+            highi: DO i=pieub(iproc),pielb(iproc),-1
+
+               DO j=MAX(1,pjelb(iproc)-jprecj), MIN(jpjglo,pjeub(iproc)+jprecj)
+
+                  IF ( depth(i,j) == 1 ) THEN
+                     ! We've found a wet point in this column so this is as far 
+                     ! as we can trim.
+                     newbound = MIN(i + jpreci + nextra, pieub(iproc))
+
+                     IF( PARTIT_DEBUG )THEN
+                        WRITE (*,"('Sub-domain ',I3,': High-i loop: col ',I3, &
+                          &   ' contains water. New bound is ',I3)") iproc,i,newbound
+                     END IF
+
+                     EXIT highi
+
+                  ENDIF
+               ENDDO
+            ENDDO highi
+
+            ! Update with the new boundary and monitor.
+
+          IF ( newbound.NE.pieub(iproc) ) THEN
+             
+             IF( TRIM_DEBUG )THEN
+                IF ( (pieub(iproc)-newbound) .GT.1 ) THEN
+                   WRITE (*,'(a,i5,3(a,i6))') ' Sub-domain ',iproc &
+                                    ,' trimmed ',pieub(iproc)-newbound &
+                          ,' cols: ',newbound+1,' to ',pieub(iproc)
+                ELSE
+                   WRITE (*,'(a,i5,3(a,i6))') ' Sub-domain ',iproc &
+                                    ,' trimmed ',pieub(iproc)-newbound &
+                                    ,' col : ',pieub(iproc)
+                ENDIF
+             END IF
+
+             pieub(iproc) = newbound
+             isTrimmed(2,iproc) = .TRUE.
+          ENDIF
+
+          ! Look at the low-j rows (Southern edge of domain)
+
+          newbound = pjelb(iproc)
+          lowj: DO j=pjelb(iproc),pjeub(iproc),1
+             
+             ! Scan this row for wet points
+
+             DO i=MAX(1,pielb(iproc)-jpreci),MIN(jpiglo,pieub(iproc)+jpreci)
+                IF ( depth(i,j) == 1) THEN
+                   newbound = MAX(j - jpreci - nextra, pjelb(iproc))
+
+                   IF( PARTIT_DEBUG )THEN
+                      WRITE (*,*) 'Sub-domain ',iproc,': Low-j loop: column ',j &
+                            ,' is land. New bound is ',newbound
+                   END IF
+
+                   EXIT lowj
+                ENDIF
+             ENDDO
+
+          ENDDO lowj
+
+
+          ! Update with the new boundary and monitor.
+
+          IF ( newbound.NE.pjelb(iproc) ) THEN
+
+             IF( TRIM_DEBUG )THEN
+                IF ( (newbound-pjelb(iproc)) .GT.1 ) THEN
+                   WRITE (*,'(a,i5,3(a,i6))') ' Sub-domain ',iproc &
+                                     ,' trimmed ',newbound-pjelb(iproc) &
+                                     ,' rows: ',pjelb(iproc),' to ',newbound-1
+                ELSE
+                   WRITE (*,'(a,i5,3(a,i6))') ' Sub-domain ',iproc &
+                                     ,' trimmed ',newbound-pjelb(iproc) &
+                                     ,' row : ',pjelb(iproc)
+                ENDIF
+             END IF
+
+             pjelb(iproc) = newbound
+             isTrimmed(3,iproc) = .TRUE.
+         ENDIF
+
+         ! Look at the high-j rows (Northern edge of domain)
+
+         newbound = pjeub(iproc)
+         highj: DO j=pjeub(iproc),pjelb(iproc),-1
+
+            ! Scan this row for wet points
+
+            DO i=MAX(1,pielb(iproc)-jpreci),MIN(jpiglo,pieub(iproc)+jpreci)
+               IF ( depth(i,j) == 1 ) THEN
+                  newbound = MIN(j + jpreci + nextra, pjeub(iproc))
+
+                  IF( TRIM_DEBUG )THEN
+                     WRITE(*,*) 'Sub-domain ',iproc,': High-j loop: column ',j &
+                                             ,' is land. New bound is ',newbound
+                  END IF
+
+                  EXIT highj
+               ENDIF
+            ENDDO
+         ENDDO highj
+
+          ! Update with the new boundary and monitor.
+
+          IF ( newbound.NE.pjeub(iproc) ) THEN
+
+             IF( TRIM_DEBUG )THEN
+                IF ( pjeub(iproc)-newbound.GT.1 ) THEN
+                   WRITE (*,'(a,i5,3(a,i6))') ' Sub-domain ',iproc &
+                                  ,' trimmed ',pjeub(iproc)-newbound &
+                                  ,' rows: ',newbound+1,' to ',pjeub(iproc)
+                ELSE
+                   WRITE (*,'(a,i5,3(a,i6))') ' Sub-domain ',iproc &
+                                  ,' trimmed ',pjeub(iproc)-newbound &
+                                  ,' row : ',pjeub(iproc)
+                ENDIF
+             ENDIF
+             pjeub(iproc) = newbound
+             isTrimmed(1,iproc) = .TRUE.
+          ENDIF
+
+          ! Reset the size of the sub-domain.
+
+          piesub(iproc) = pieub(iproc)-pielb(iproc)+1
+          pjesub(iproc) = pjeub(iproc)-pjelb(iproc)+1
+
+          ! endif active_subdomain
+
+       ENDDO subdomain_loop
+
+   END SUBROUTINE part_trim
+
+
+   SUBROUTINE finish_partition(imask)
+      IMPLICIT NONE
+      INTEGER, DIMENSION(:,:), INTENT(in) :: imask
+      ! Locals
+      INTEGER :: iproc, ierr
+
+      ! Set the external boundary flags before boundaries are
+      ! altered by the trimming process and it becomes more difficult
+      ! to recognize which were the external boundaries.
+      
+      DO iproc=1, nprocp, 1
+         pilbext(iproc) = pielb(iproc).EQ.1
+         piubext(iproc) = pieub(iproc).EQ.jpiglo
+         pjlbext(iproc) = pjelb(iproc).EQ.1
+         pjubext(iproc) = pjeub(iproc).EQ.jpjglo
+      ENDDO
+
+      ! Trim off redundant rows and columns containing all land.
+
+      IF ( trim_land ) THEN
+         CALL part_trim ( imask, trimmed, ierr )
+      ELSE
+         trimmed(1:4,1:nprocp) = .FALSE.
+      ENDIF
+
+!!$
+!!$      ! Lower boundary (long.) of sub-domain, GLOBAL coords
+!!$      ! before correction for global halos
+!!$      nimpp = pielb(narea)
+!!$
+!!$      ! Is the domain on an external LONGITUDE boundary?
+!!$      nbondi = 0
+!!$      ilbext = pilbext(narea)
+!!$      IF(ilbext)THEN
+!!$         nbondi = -1
+!!$      END IF
+!!$
+!!$      IF( (.NOT. ilbext) .OR. (ilbext .AND. trimmed(widx,narea)) )THEN
+!!$         ! It isn't, which means we must shift its lower boundary by
+!!$         ! -jpreci to allow for the overlap of this domain with its
+!!$         ! westerly neighbour.
+!!$         nimpp = nimpp - jpreci
+!!$      END IF
+!!$
+!!$      iubext = piubext(narea)
+!!$      IF(iubext)THEN
+!!$         nbondi = 1
+!!$         IF(ilbext)nbondi = 2 ! Both East and West boundaries are at 
+!!$                              ! edges of global domain
+!!$      END IF
+!!$
+!!$      ! Set local values for limits in global coords of the sub-domain 
+!!$      ! owned by this process.
+!!$      ielb   = pielb (narea)
+!!$      ieub   = pieub (narea)
+!!$      iesub  = piesub(narea)
+!!$
+!!$      jpi  = iesub + 2*jpreci ! jpi is the same for all domains - this is
+!!$                              ! what original decomposition did
+!!$      nlci = jpi
+!!$
+!!$      ! If the domain is at the edge of the model domain and a cyclic 
+!!$      ! East-West b.c. is in effect then it already incorporates one 
+!!$      ! extra halo column (because of the way the model domain itself is
+!!$      ! set-up). If we've trimmed-off dry rows and columns then, even if
+!!$      ! a domain is on the model boundary, it may still need a halo so
+!!$      ! we add one.
+!!$      IF( nbondi == -1 .AND. (.NOT. trimmed(widx,narea))  )THEN
+!!$           ! Western boundary
+!!$           ! First column of global domain is actually a halo but NEMO
+!!$           ! still sets nldi to 1.
+!!$           nldi = 1            ! Lower bnd of int. region of sub-domain, LOCAL
+!!$           nlei = nldi + iesub - 1
+!!$           nlci = nlei + jpreci
+!!$           jpi  = nlci
+!!$
+!!$        ELSE IF( nbondi == 1 .AND. (.NOT. trimmed(eidx,narea)) )THEN
+!!$           ! Eastern boundary
+!!$           nldi = 1 + jpreci
+!!$           ! Last column of global domain is actually a halo
+!!$           nlei = nldi + iesub - 1
+!!$           nlci = nlei
+!!$           jpi  = nlei
+!!$
+!!$        ELSE IF( nbondi == 2)THEN
+!!$           ! Both boundaries are on the edges of the global domain
+!!$           IF(trimmed(widx,narea))THEN
+!!$              nldi = 1 + jpreci
+!!$           ELSE
+!!$              nldi = 1
+!!$           END IF
+!!$           nlei = nldi + iesub - 1
+!!$
+!!$           IF(trimmed(eidx,narea))THEN
+!!$              nlci = nlei + jpreci
+!!$           ELSE
+!!$              nlci = nlei
+!!$           END IF
+!!$           jpi  = nlci
+!!$
+!!$        ELSE
+!!$           ! We have no external boundaries to worry about
+!!$           nldi = 1 + jpreci 
+!!$           nlei = nldi + iesub - 1 !
+!!$        END IF
+!!$
+!!$        jpim1 = jpi - 1
+!!$
+!!$        ! Lower ext. boundary (lat.) of sub-domain, global coords
+!!$        njmpp= pjelb (narea)
+!!$
+!!$        ! Is the domain on an external LATITUDE boundary?
+!!$        nbondj = 0
+!!$        jlbext = pjlbext(narea)
+!!$        IF(jlbext)THEN
+!!$           nbondj = -1
+!!$        ENDIF
+!!$
+!!$        IF((.NOT. jlbext) .OR. (jlbext .AND. trimmed(sidx,narea)) )THEN
+!!$           ! It isn't, which means we must shift its lower boundary by
+!!$           ! -jprecj to allow for the overlap of this domain with its
+!!$           ! southerly neighbour.
+!!$           njmpp = njmpp - jprecj
+!!$        END IF
+!!$        jubext = pjubext(narea)
+!!$        IF(jubext)THEN
+!!$           nbondj = 1
+!!$           IF(jlbext)nbondj = 2
+!!$        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
+!!$
+!!$      ! Unlike the East-West boundaries, the global domain does not include
+!!$      ! halo (rows) at the Northern and Southern edges. In fact, there is no
+!!$      ! cyclic boundary condition in the North-South direction so there are no
+!!$      ! halos at all at the edges of the global domain.
+!!$      IF( nbondj == -1 .AND. (.NOT. trimmed(sidx,narea)) )THEN
+!!$         ! Southern edge
+!!$         nldj = 1                ! Start of internal region (local coords)
+!!$         nlej = nldj + jesub - 1 ! Upper bound of int. region of sub-domain, local
+!!$         nlcj = nlej + jprecj
+!!$         jpj  = nlcj
+!!$      ELSE IF( nbondj ==  1 .AND. (.NOT. trimmed(nidx,narea)) )THEN
+!!$         ! Northern edge
+!!$         nldj = 1+jprecj       ! Start of internal region (local coords)
+!!$         nlej = nldj + jesub - 1
+!!$         nlcj = nlej
+!!$         jpj  = nlcj
+!!$         jpj = jpj + 1 ! Add one extra row for zero BC along N edge as
+!!$                       ! happens in orig. code when MOD(jpjglo,2)!=0
+!!$                       ! Many loops go up to j=jpjm1 so unless jpj>nlej
+!!$                       ! they won't cover the whole domain.
+!!$      ELSE IF( nbondj == 2)THEN
+!!$         ! Both boundaries are on the edges of the global domain
+!!$
+!!$         IF( trimmed(sidx,narea) )THEN
+!!$            nldj = 1+jprecj
+!!$         ELSE
+!!$            nldj = 1
+!!$         END IF
+!!$         nlej = nldj + jesub - 1
+!!$
+!!$         IF( trimmed(nidx,narea) )THEN
+!!$            nlcj = nlej + jprecj
+!!$            jpj  = nlcj
+!!$         ELSE
+!!$            nlcj = nlej
+!!$            jpj  = nlcj
+!!$         END IF
+!!$         jpj = jpj + 1 ! Add one extra row for zero BC along N edge as
+!!$                       ! happens in orig. code when MOD(jpjglo,2)!=0
+!!$      ELSE
+!!$         ! We have no external boundaries to worry about
+!!$         nldj = 1+jprecj    ! Lower bound of int. region of sub-domain, local
+!!$         nlej = nldj + jesub - 1
+!!$      END IF
+!!$
+!!$      jpjm1 = jpj-1
+!!$      jpij  = jpi*jpj 
+
+   END SUBROUTINE finish_partition
+
+
+   SUBROUTINE eval_partition( cores_per_node, 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)
+      !   + ...
+      ! }
+      IMPLICIT NONE
+      !     Arguments
+      INTEGER, INTENT(in) :: cores_per_node
+      INTEGER, INTENT(in) :: nx, ny
+      INTEGER, INTENT(in) :: mask(nx,ny)
+      REAL(wp), DIMENSION(pv_num_scores), INTENT(out) :: score
+      !     Local variables
+      INTEGER :: iproc, inode, i, j
+
+      REAL(wp) :: score_wet, score_dry
+      REAL(wp) :: score_pli, score_plx
+      REAL(wp) :: score_pci, score_pcx
+      REAL(wp) :: score_nli, score_nlx
+      REAL(wp) :: score_nci, score_ncx
+      REAL(wp) :: score_tli, score_tlx
+      REAL(wp) :: score_tci, score_tcx
+      
+      REAL(wp) :: score_too_narrow
+      
+      REAL(wp) :: proc_overall_score
+      REAL(wp) :: proc_comm_score
+      REAL(wp) :: node_comm_score
+
+      REAL(wp), PARAMETER :: weight_wet  =  1.00D0
+      REAL(wp), PARAMETER :: weight_dry  =  0.9D0
+      
+      REAL(wp), PARAMETER :: weight_pli  =  0.01D0
+      REAL(wp), PARAMETER :: weight_plx  =  0.20D0
+      REAL(wp), PARAMETER :: weight_pci  =  0.50D0
+      REAL(wp), PARAMETER :: weight_pcx  = 10.00D0
+      
+      REAL(wp), PARAMETER :: weight_nli  =  0.00D0
+      REAL(wp), PARAMETER :: weight_nlx  =  0.20D0
+      REAL(wp), PARAMETER :: weight_nci  =  0.00D0
+      REAL(wp), PARAMETER :: weight_ncx  = 10.00D0
+        
+      REAL(wp), PARAMETER :: weight_tli  =  0.00D0
+      REAL(wp), PARAMETER :: weight_tlx  =  0.01D0
+      REAL(wp), PARAMETER :: weight_tci  =  0.00D0
+      REAL(wp), PARAMETER :: weight_tcx  =  0.50D0
+
+      INTEGER :: peer, last_peer
+
+      ! Which node is each process on?
+      ! Assumes that first cores_per_node ranks are assigned to node 0,
+      ! next cores_per_node ranks are assigned to node 1, etc
+      INTEGER, ALLOCATABLE :: node(:)
+
+      ALLOCATE(node(nprocp))
+
+      DO iproc=1, nprocp
+         node(iproc) = (iproc-1)/cores_per_node
+      END DO
+
+      ! Calculate maximum per processor score
+
+      score(:) = 0.0_wp
+
+      ! Total (over all processors) off node comms 
+      score_tli  = 0.0_wp
+      score_tlx  = 0.0_wp
+      score_tci  = 0.0_wp
+      score_tcx  = 0.0_wp
+      
+      ! Set this to pv_awful if any sub-domain is too narrow.
+      score_too_narrow = 0.0_wp
+
+      ! loop over nodes in job, counting from 0
+      node_loop: DO inode=0, (nprocp-1)/cores_per_node
+
+        score_nli  = 0.0_wp
+        score_nlx  = 0.0_wp
+        score_nci  = 0.0_wp
+        score_ncx  = 0.0_wp
+
+        ! loop over processes in the node
+        proc_loop: DO iproc=1+inode*cores_per_node, &
+                            MIN(nprocp,(inode+1)*cores_per_node) 
+
+           score_wet  = DBLE(pnactive(iproc))
+           score_dry  = DBLE(piesub(iproc)*pjesub(iproc)-score_wet)
+
+           score_pli  = 0.0_wp
+           score_plx  = 0.0_wp
+           score_pci  = 0.0_wp
+           score_pcx  = 0.0_wp
+
+           ! Things sometimes go wrong when a sub-domain has very
+           ! narrow partitions (2 grid points or less).
+           ! Prevent such problematic partitions from being selected.
+           IF (      ((pieub(iproc)-pielb(iproc)) .LE. 2) &
+                .OR. ((pjeub(iproc)-pjelb(iproc)) .LE. 2) ) THEN
+              score_too_narrow = pv_awful
+           END IF
+
+           IF (.NOT. pjlbext(iproc)) THEN
+              j=pjelb(iproc)
+              IF (j .GT. 1) THEN
+                 last_peer = -1
+                 DO i=pielb(iproc),pieub(iproc)
+                    IF (       (mask(i,j) .NE. land) &
+                         .AND. (mask(i,j-1) .NE. land)) THEN
+                       peer=iprocmap(i,j-1)
+                       ! Total points involved in messages
+                       IF (node(peer) .EQ. inode) THEN
+                          score_pli = score_pli+1.0_wp
+                       ELSE
+                          score_plx = score_plx+1.0_wp
+                       END IF
+                       ! Total number of messages
+                       IF (peer .NE. last_peer) THEN
+                          last_peer = peer
+                          IF (node(peer) .EQ. inode) THEN
+                             score_pci = score_pci+1.0_wp
+                          ELSE
+                             score_pcx = score_pcx+1.0_wp
+                          END IF
+                       END IF
+                    END IF
+                 END DO
+              END IF
+           END IF
+  
+           IF (.NOT. pjubext(iproc)) THEN
+              j=pjeub(iproc)
+              IF (j .LT. ny) THEN
+                 last_peer = -1
+                 DO i=pielb(iproc),pieub(iproc)
+                    IF (      (mask(i,j)   .NE. land)  &
+                        .AND. (mask(i,j+1) .NE. land)) THEN
+                       peer=iprocmap(i,j+1)
+                       ! Total points involved in messages
+                       IF (node(peer) .EQ. inode) THEN
+                          score_pli = score_pli+1.0_wp
+                       ELSE
+                          score_plx = score_plx+1.0_wp
+                       END IF
+                       ! Total number of messages
+                       IF (peer .NE. last_peer) THEN
+                          last_peer = peer
+                          IF (node(peer) .EQ. inode) THEN
+                             score_pci = score_pci+1.0_wp
+                          ELSE
+                             score_pcx = score_pcx+1.0_wp
+                          END IF
+                       END IF
+                    END IF
+                 END DO
+              END IF
+           END IF
+           
+           IF (.NOT. pilbext(iproc)) THEN
+              i=pielb(iproc)
+              IF (i .GT. 1) THEN
+                 last_peer = -1
+                 DO j=pjelb(iproc),pjeub(iproc)
+                    IF (      (mask(i,j)   .NE. land) &
+                        .AND. (mask(i-1,j) .NE. land)) THEN
+                       peer=iprocmap(i-1,j)
+                       ! Total points involved in messages
+                       IF (node(peer) .EQ. inode) THEN
+                          score_pli = score_pli+1.0_wp
+                       ELSE
+                          score_plx = score_plx+1.0_wp
+                       END IF
+                       ! Total number of messages
+                       IF (peer .NE. last_peer) THEN
+                          last_peer = peer
+                          IF (node(peer) .EQ. inode) THEN
+                             score_pci = score_pci+1.0_wp
+                          ELSE
+                             score_pcx = score_pcx+1.0_wp
+                          END IF
+                       END IF
+                    END IF
+                 END DO
+              END IF
+           END IF
+
+          IF (.NOT. piubext(iproc)) THEN
+            i=pieub(iproc)
+            IF (i .LT. nx) THEN
+              last_peer = -1
+              DO j=pjelb(iproc),pjeub(iproc)
+                IF (      (mask(i,j)   .NE. land)  &
+                    .AND. (mask(i+1,j) .NE. land)) THEN
+                  peer=iprocmap(i+1,j)
+                  ! Total points involved in messages
+                  IF (node(peer) .EQ. inode) THEN
+                    score_pli = score_pli+1.0_wp
+                  ELSE
+                    score_plx = score_plx+1.0_wp
+                  END IF
+                  ! Total number of messages
+                  IF (peer .NE. last_peer) THEN
+                    last_peer = peer
+                    IF (node(peer) .EQ. inode) THEN
+                      score_pci = score_pci+1.0_wp
+                    ELSE
+                      score_pcx = score_pcx+1.0_wp
+                    END IF
+                  END IF
+                END IF
+              END DO
+            END IF
+          END IF
+
+          score_nli = score_nli + score_pli
+          score_nlx = score_nlx + score_plx
+          score_nci = score_nci + score_pci
+          score_ncx = score_ncx + score_pcx
+
+          proc_overall_score = weight_wet*score_wet &
+                            + weight_dry*score_dry  &
+                            + weight_pli*score_pli  &
+                            + weight_plx*score_plx  &
+                            + weight_pci*score_pci  &
+                            + weight_pcx*score_pcx
+
+          proc_comm_score    = weight_pli*score_pli &
+                             + weight_plx*score_plx &
+                             + weight_pci*score_pci &
+                             + weight_pcx*score_pcx
+
+          IF (score(pv_index_overall) < proc_overall_score) &
+              score(pv_index_overall) = proc_overall_score
+          IF (score(pv_index_pcomms) < proc_comm_score)     &
+              score(pv_index_pcomms) = proc_comm_score
+          IF (score(pv_index_wet) < score_wet) &
+              score(pv_index_wet) = score_wet
+          IF (score(pv_index_dry) < score_dry) &
+              score(pv_index_dry) = score_dry
+          IF (score(pv_index_pli) < score_pli) &
+              score(pv_index_pli) = score_pli
+          IF (score(pv_index_plx) < score_plx) &
+              score(pv_index_plx) = score_plx
+          IF (score(pv_index_pci) < score_pci) &
+              score(pv_index_pci) = score_pci
+          IF (score(pv_index_pcx) < score_pcx) &
+              score(pv_index_pcx) = score_pcx
+  
+        END DO proc_loop
+
+        score_tli = score_tli + score_nli
+        score_tlx = score_tlx + score_nlx
+        score_tci = score_tci + score_nci
+        score_tcx = score_tcx + score_ncx
+
+        node_comm_score    = weight_nli*score_nli &
+                           + weight_nlx*score_nlx &
+                           + weight_nci*score_nci &
+                           + weight_ncx*score_ncx
+
+        IF (score(pv_index_ncomms) < node_comm_score) &
+            score(pv_index_ncomms) = node_comm_score
+        IF (score(pv_index_nli) < score_nli) &
+            score(pv_index_nli) = score_nli
+        IF (score(pv_index_nlx) < score_nlx) &
+            score(pv_index_nlx) = score_nlx
+        IF (score(pv_index_nci) < score_nci) &
+            score(pv_index_nci) = score_nci
+        IF (score(pv_index_ncx) < score_ncx) &
+            score(pv_index_ncx) = score_ncx
+
+      END DO node_loop
+
+      score(pv_index_tcomms) = weight_tli*score_tli &
+                             + weight_tlx*score_tlx &
+                             + weight_tci*score_tci &
+                             + weight_tcx*score_tcx
+      
+      score(pv_index_tli) = score_tli
+      score(pv_index_tlx) = score_tlx
+      score(pv_index_tci) = score_tci
+      score(pv_index_tcx) = score_tcx
+
+      score(pv_index_overall) = score(pv_index_overall) &
+                              + score(pv_index_ncomms)  &
+                              + score(pv_index_tcomms)  &
+                              + score_too_narrow
+
+      DEALLOCATE(node)
+
+    END SUBROUTINE eval_partition
+
+
+    SUBROUTINE calc_perms( nfactors, factors,     &
+                           ndf, df, multiplicity, &
+                           nperms )
+      IMPLICIT NONE
+
+      !     Subroutine arguments
+      ! Number of factors (including repetitions)
+      INTEGER, INTENT(in) :: nfactors
+      ! Factors (in descending order)
+      INTEGER, INTENT(in) :: factors(nfactors)
+      ! Number of distinct factors
+      INTEGER, INTENT(out) :: ndf
+      ! Distinct factors (in ascending order)
+      INTEGER :: df(nfactors)
+      ! Multiplicity of each distinct factor
+      INTEGER :: multiplicity(nfactors)
+      ! Number of distinct permutations
+      INTEGER, INTENT(out) :: nperms
+
+      ! Local variables
+      INTEGER :: product
+      INTEGER :: i, j
+
+      product = factors(nfactors)
+      ndf = 1
+      df(:) = 0
+      df(1) = factors(nfactors)
+      multiplicity(:) = 0
+      multiplicity(1) = 1
+      
+      DO i=nfactors-1,1,-1
+        product = product*factors(i)
+        IF (factors(i) .EQ. df(ndf)) THEN
+          multiplicity(ndf) = multiplicity(ndf)+1
+        ELSE
+          ndf = ndf+1
+          df(ndf) = factors(i)
+          multiplicity(ndf) = 1
+        END IF
+      END DO
+
+!     A careful code would try to avoid overflow here
+      nperms = 1
+      DO i=1, nfactors
+        nperms = nperms*i
+      END DO
+      DO i=1, ndf
+        DO j=1, multiplicity(i)
+          nperms = nperms/j
+        END DO
+      END DO
+
+!     At this point, nperms is the number of distinct permutations
+!     of the factors provided. But each of these permutations can
+!     be split between x and y in (nfactors+1) ways, e.g.
+!       x=(2,2,3), y=()
+!       x=(2,3),   y=(2)
+!       x=(3),     y=(2,2)
+!       x=(),      y=(2,2,3)
+
+      nperms = nperms*(nfactors+1)
+
+      WRITE (*,FMT="('Distinct factorisations: ', I4/)") nperms
+          
+    END SUBROUTINE calc_perms
+
+
+    
+    SUBROUTINE get_perm_factors( iperm, nf, ndf, df, m, &
+                                 fx, nfx, fy, nfy,      &
+                                 prodx, prody )
+         IMPLICIT NONE
+         !!------------------------------------------------------------------
+         !     Our goal is to visit a particular permutation, 
+         !     number perm ( 0 <= perm <= nperms-1 )
+         !     of nf things, of ndf distinct values (actually the prime
+         !     factors of number of MPI tasks), each of which can be repeated
+         !     with multiplicity m_i
+         !     assert nf = sum_{i=1..ndf} m(i)
+         !     
+         !     We don't care about lexicographic ordering, but we do
+         !     need to ensure that we don't visit any permutation twice,
+         !     in a loop over 0..nperms-1.
+         !     Textbook methods typically assume that all the things being
+         !     permuted are distinct.
+         !
+         !     We use what I call a nested variable radix method.
+         !
+         !     Stephen Pickles, STFC
+         !     Taken from POLCOMS code by Andrew Porter.
+         !!------------------------------------------------------------------
+         !     Arguments
+         INTEGER, INTENT(in)  :: iperm, nf, ndf
+         INTEGER, INTENT(in), DIMENSION(ndf) :: df, m
+         INTEGER, INTENT(out), DIMENSION(nf) :: fx, fy
+         INTEGER, INTENT(out) :: nfx, nfy
+         INTEGER, INTENT(out) :: prodx, prody
+         !
+         !     Local variables
+         !
+         INTEGER :: perm, split
+         INTEGER, DIMENSION(nf) :: bin, a
+         INTEGER, DIMENSION(ndf) :: ways
+         INTEGER, DIMENSION(0:ndf) :: root, representation
+         INTEGER :: i, j, k, v, p, q, r
+         INTEGER :: unfilled, pm
+         INTEGER, SAVE :: myinst = 0
+         LOGICAL, PARAMETER :: debug=.FALSE.
+         !!------------------------------------------------------------------
+
+         ! This is a serial version so instead of using MPI rank
+         ! use the count of how many times we've been called
+         myinst = myinst + 1
+
+         perm = iperm / (nf+1)
+         !     Where to split between x and y 
+         split = MOD(iperm, (nf+1))
+
+         !     interpret ways(i) is the number of ways of distributing 
+         !     m(i) copies of the i'th distinct factor into the remaining
+         !     bins
+         k = nf
+         DO i=1,ndf
+            ways(i) = k
+            DO j=2,m(i)
+               ways(i) = ways(i)*(k-j+1)/j
+            END DO
+            k = k-m(i)
+         END DO
+
+         !     compute (outer) radices
+         !     root is the variable radix basis corresponding to ways
+         !     root(ndf) is always 1
+         root(ndf) = 1
+         root(0:ndf-1) = ways(1:ndf)
+         DO i=ndf-1,0,-1
+            root(i)=root(i)*root(i+1)
+         END DO
+
+         bin(:) = 0
+         unfilled = nf
+
+         r = perm
+         !     Next line is valid as long as perm < nperms
+         representation(0) = 0
+         DO i=1, ndf
+            p = r/root(i)
+            r = r - p*root(i)
+            representation(i) = p
+
+            !       At this point, we are considering distinct ways to
+            !       distribute m(i) copies of the i'th distinct factor into
+            !       the unfilled remaining bins. We want to select the p'th one.
+
+            DO j=1,unfilled
+               a(j) = j
+            END DO
+            q = 0
+            find_p: DO
+               IF (q .GE. p) EXIT find_p
+
+               k=m(i)
+               k_loop: DO
+                  IF (a(k) .EQ. (unfilled - m(i) + k)) THEN
+                     k=k-1
+                  ELSE
+                     EXIT k_loop
+                  END IF
+               END DO k_loop
+               a(k) = a(k) + 1
+               DO v=k+1,m(i)
+                  a(v) = a(k) + v - k
+               END DO
+               q = q+1
+            END DO find_p
+
+            IF (debug) THEN
+               WRITE (*,'(A3)',advance='no') 'a=('
+               DO j=1, m(i)-1
+                  WRITE (*,'(I3,A1)',advance='no') a(j), ','
+               END DO
+               WRITE (*,'(I3,A1)') a(m(i)), ')'
+            END IF
+
+            DO j=1, m(i)
+               pm=a(j)-j+1
+
+               ! put this factor in the pm'th empty bin
+               v = 1
+               find_bin: DO k=1, nf
+                  IF (bin(k) .EQ. 0) THEN
+                     IF (v .EQ. pm) THEN
+                        bin(k) = df(i)
+                        unfilled = unfilled-1
+                        EXIT find_bin
+                     ELSE
+                        v=v+1
+                     END IF
+                  END IF
+               END DO find_bin
+
+            END DO
+            
+         END DO
+
+         !     We have identified the (perm)th distinct permutation,
+         !     but we still need to split the factors between x and y.
+         fx(:) = 0
+         prodx = 1
+         DO i=1,split
+            fx(i)=bin(i)
+            prodx=prodx*fx(i)
+         END DO
+         nfx=split
+
+         fy(:) = 0 
+         prody = 1
+         j=1
+         DO i=split+1,nf
+            fy(j)=bin(i)
+            prody=prody*fy(j)
+            j=j+1
+         END DO
+         nfy=nf-nfx
+
+            WRITE (6,'(A14,I6,A1,I6,A2)',advance='no') &
+                      'factorisation[', myinst, ']', iperm, ' ('
+            DO k=1,ndf-1
+               WRITE (6,'(I4,A1)',advance='no') representation(k), ','
+            END DO
+            WRITE (6,'(I4,A1)',advance='no') representation(ndf), ')'
+     
+            CALL print_factors(nfx,fx,nfy,fy)
+
+      END SUBROUTINE get_perm_factors
+
+
+      SUBROUTINE print_factors(nfx,fx,nfy,fy)
+         IMPLICIT NONE
+         INTEGER, INTENT(in) :: nfx, nfy
+         INTEGER, INTENT(in) :: fx(nfx), fy(nfy)
+         INTEGER :: k
+
+         IF (nfx .GT. 0) THEN
+            WRITE (*,'(A1)',advance='no') ' '
+            DO k=1,nfx-1
+               WRITE (*,'(I4,A1)',advance='no') fx(k), ','
+            END DO
+            WRITE (*,'(I4)',advance='no') fx(nfx)
+         END IF
+         WRITE (*,'(A1)',advance='no') ':'
+         IF (nfy .GT. 0) THEN
+            DO k=1,nfy-1
+               WRITE (*,'(I4,A1)',advance='no') fy(k), ','
+            END DO
+            WRITE (*,'(I4)',advance='no') fy(nfy)
+         END IF
+         WRITE (*,'(A1)') ' '
+
+      END SUBROUTINE print_factors
+
+
+      CHARACTER(len=16) FUNCTION num_to_string(number)
+         INTEGER, INTENT(in) :: number
+
+         CHARACTER*16 :: outs
+      
+         WRITE (outs,'(I15)') number
+         num_to_string = ADJUSTL(outs)
+
+      END FUNCTION num_to_string
+
+
+      SUBROUTINE factor_string(fstr,nfx,fx,nfy,fy)
+         IMPLICIT NONE
+         CHARACTER*256, INTENT(out) :: fstr
+         INTEGER, INTENT(in) :: nfx, nfy
+         INTEGER, INTENT(in) :: fx(nfx), fy(nfy)
+         ! Locals
+         INTEGER :: k
+
+         fstr = ' '
+         IF (nfx .GT. 0) THEN
+            DO k=1,nfx-1
+               fstr = TRIM(fstr)//TRIM(num_to_string(fx(k)))//'x'
+            END DO
+            fstr = TRIM(fstr)//TRIM(num_to_string(fx(nfx)))
+         END IF
+         fstr = TRIM(fstr)//'-'
+         IF (nfy .GT. 0) THEN
+            DO k=1,nfy-1
+               fstr = TRIM(fstr)//TRIM(num_to_string(fy(k)))//'x'
+            END DO
+            fstr = TRIM(fstr)//TRIM(num_to_string(fy(nfy)))
+         END IF
+      END SUBROUTINE factor_string
+
+      FUNCTION iprocmap(ia, ja)
+!!------------------------------------------------------------------
+! Returns the process number (1...N) of the process whose sub-domain
+! contains the point with global coordinates (i,j).
+! If no process owns the point, returns zero.
+
+!     i                       int   input     global x-coordinate
+!     j                       int   input     global y-coordinate
+
+!         Mike Ashworth, CLRC Daresbury Laboratory, July 1999
+!         Andrew Porter, STFC Daresbury Laboratory, May  2008
+!!------------------------------------------------------------------
+        IMPLICIT NONE
+
+        !  Function arguments.
+        INTEGER             :: iprocmap
+        INTEGER, INTENT(in) :: ia, ja
+        ! Local variables.
+        INTEGER :: iproc, i, j, iwidth
+
+        iprocmap = 0
+
+        ! Make sure we don't change variable values in calling
+        ! routine...
+        i = ia
+        j = ja
+        IF(cyclic_bc)THEN
+           ! Allow for fact that first and last columns in global domain
+           ! are actually halos
+           iwidth = jpiglo - 2*jpreci
+           IF(i >= jpiglo) i = ia - iwidth
+           IF(i <= 1     ) i = ia + iwidth
+           ! No cyclic BCs in North/South direction
+           !IF(j > jpjglo) j = ja - jpjglo
+           !IF(j < 1     ) j = ja + jpjglo
+        END IF
+
+        ! Search the processes for the one owning point (i,j).
+
+        DO iproc=1,nprocp
+           IF ( pielb(iproc).LE.i .AND. i.LE.pieub(iproc) .AND. &
+                pjelb(iproc).LE.j .AND. j.LE.pjeub(iproc) ) THEN
+              iprocmap = iproc
+              EXIT
+           ENDIF
+        ENDDO
+     END FUNCTION iprocmap
+
+    !===================================================================
+
+    SUBROUTINE write_partition(imask)
+      IMPLICIT none
+      INTEGER, DIMENSION(:,:), INTENT(in) :: imask
+      ! Locals
+      INTEGER, PARAMETER :: funit = 1099
+      INTEGER :: i, j, iproc, ierr
+      INTEGER :: idom
+      CHARACTER(LEN=4) :: intStr
+      ! Space between plot and (each) axis
+      INTEGER, PARAMETER :: axis_offset = 4
+      INTEGER, PARAMETER :: tic_len = 4
+      INTEGER, PARAMETER :: tic_label_space = 2
+      INTEGER, PARAMETER :: xtic_incr = 50
+      INTEGER, PARAMETER :: ytic_incr = 50
+
+      ! Write the file that NEMO can read
+
+      OPEN(UNIT=funit, file='partition.dat.new', status='REPLACE', &
+           ACTION='WRITE', IOSTAT=ierr)
+      IF(ierr /= 0)THEN
+         WRITE(*,*) 'ERROR: write_partition: failed to write partition description to file.'
+         RETURN
+      END IF
+
+      WRITE(UNIT=funit,FMT="(I5,1x,I5)") jpni, jpnj
+
+      DO idom = 1,nprocp, 1
+         WRITE(UNIT=funit,FMT="('Domain: ',I5)") idom
+         WRITE(UNIT=funit,FMT="(I5,1x,I5)") pielb(idom), pjelb(idom)
+         WRITE(UNIT=funit,FMT="(I5,1x,I5)") pieub(idom), pjeub(idom)
+         WRITE(UNIT=funit,FMT="(L5,1x,L5)") pilbext(idom), trimmed(widx,idom)
+         WRITE(UNIT=funit,FMT="(L5,1x,L5)") piubext(idom), trimmed(eidx,idom)
+         WRITE(UNIT=funit,FMT="(L5,1x,L5)") pjlbext(idom), trimmed(sidx,idom)
+         WRITE(UNIT=funit,FMT="(L5,1x,L5)") pjubext(idom), trimmed(nidx,idom)
+      END DO
+
+      CLOSE(UNIT=funit)
+
+      ! Write out domains in postscript
+
+      OPEN(UNIT=997, FILE="domain_decomp.ps", &
+           STATUS='REPLACE', ACTION='WRITE', IOSTAT=ierr)
+
+      IF(ierr /= 0)THEN
+         WRITE(*,*) 'WARNING: failed to open domain_decomp.ps for writing'
+         RETURN
+      END IF
+
+      ! Header of ps file
+      WRITE (997,FMT='("%!PS-Adobe-1.0")')
+      WRITE (997,FMT='("/Helvetica findfont %load the font dictionary")')
+      WRITE (997,FMT='("12 scalefont        %scale to 12pt")')
+      WRITE (997,FMT='("setfont             %make this the current font")')
+      WRITE (997,FMT='("/u { 2 mul } def    %set up a scaling factor")')
+
+      ! Put green circles at dry points
+      WRITE (997,FMT='("% Filled green circles at dry points")')
+      WRITE (997,FMT='("0.1 setlinewidth")') ! Thin green line
+      WRITE (997,FMT='("0 1 0 setrgbcolor")')
+      WRITE (997,FMT='("newpath")')
+      DO j = 1,jpjglo,1
+         DO i = 1,jpiglo,1
+            IF(imask(i,j) /= 1)THEN
+               ! Radius of each circle is unity
+               WRITE (997,FMT='(I3," u ",I3," u 1  0 360 arc closepath")') i, j
+               ! Use 'fill' instead of 'stroke' to get filled circles
+               WRITE (997,FMT='("fill")')
+            END IF
+         END DO
+      END DO
+
+      ! Draw the outline of the global domain in red
+      WRITE (997,FMT='("% Draw the outline of the global domain in red")')
+      WRITE (997,FMT='("3.0 setlinewidth")') ! Fat line of 3mm for global dom.
+      WRITE (997,FMT='("1 0 0 setrgbcolor")')
+      WRITE (997,FMT='("newpath")')
+      WRITE (997,FMT='("% Cursor initialization")')
+      WRITE (997,FMT='(I3," u ",1x,I3," u moveto")') 1,1
+      WRITE (997,FMT='("% Drawing the rectangle")')
+      WRITE (997,FMT='(I3," u ",1x,I3," u lineto")') jpiglo,1
+      WRITE (997,FMT='(I3," u ",1x,I3," u lineto")') jpiglo,jpjglo
+      WRITE (997,FMT='(I3," u ",1x,I3," u lineto")') 1,jpjglo
+      WRITE (997,FMT='("closepath stroke")')
+
+      ! Axes 
+      WRITE (997,FMT='("0 0 0 setrgbcolor")')! back to black
+      WRITE (997,FMT='("% Draw axes")')
+      WRITE (997,FMT='("newpath")')
+      WRITE (997,FMT='("% y axis")')
+      WRITE (997,FMT='(I3," u ",1x,I3," u moveto")') jpiglo+axis_offset,1
+      WRITE (997,FMT='(I3," u ",1x,I3," u lineto")') jpiglo+axis_offset,jpjglo
+      WRITE (997,FMT='("closepath stroke")')
+      ! Tic marks
+      ! At y==1
+      j = 1
+      WRITE (997,FMT='(I3," u ",1x,I3," u moveto")') jpiglo+axis_offset,j
+      WRITE (997,FMT='(I3," u ",1x,I3," u lineto")') jpiglo+axis_offset+tic_len,j
+      WRITE (997,FMT='("closepath stroke")')
+      ! Tic label
+      WRITE (997,FMT='(I3," u ",1x,I3," u moveto")') jpiglo+axis_offset+tic_len+tic_label_space,j
+      WRITE (997,FMT='("(1) dup stringwidth pop 2 div neg 0 rmoveto show")')
+
+      DO j = ytic_incr, jpjglo, ytic_incr
+         WRITE (997,FMT='(I3," u ",1x,I3," u moveto")') jpiglo+axis_offset,j
+         WRITE (997,FMT='(I3," u ",1x,I3," u lineto")') jpiglo+axis_offset+tic_len,j
+         WRITE (997,FMT='("closepath stroke")')
+         ! Tic label
+         WRITE (997,FMT='(I3," u ",1x,I3," u moveto")') jpiglo+axis_offset+tic_len+tic_label_space,j
+         WRITE (intStr, FMT='(I4)') j
+         WRITE (997,FMT='("(",(A),") show")') TRIM(ADJUSTL(intStr))
+       END DO
+
+      WRITE (997,FMT='("% x axis")')
+      WRITE (997,FMT='(I3," u ",1x,I3," u moveto")') 1,jpjglo+axis_offset
+      WRITE (997,FMT='(I3," u ",1x,I3," u lineto")') jpiglo,jpjglo+axis_offset
+      WRITE (997,FMT='("closepath stroke")')
+      WRITE (997,FMT='("closepath stroke")')
+      ! Tic marks
+      ! For x==1
+      i = 1
+      WRITE (997,FMT='(I3," u ",1x,I3," u moveto")') i,jpjglo+axis_offset
+      WRITE (997,FMT='(I3," u ",1x,I3," u lineto")') i,jpjglo+axis_offset+tic_len
+      WRITE (997,FMT='("closepath stroke")')
+      ! Tic label
+      WRITE (997,FMT='(I3," u ",1x,I3," u moveto")') i,jpjglo+axis_offset+tic_len+tic_label_space
+      WRITE (997,FMT='("(1) dup stringwidth pop 2 div neg 0 rmoveto show")')
+
+      DO i = xtic_incr, jpiglo, xtic_incr
+         WRITE (997,FMT='(I3," u ",1x,I3," u moveto")') i,jpjglo+axis_offset
+         WRITE (997,FMT='(I3," u ",1x,I3," u lineto")') i,jpjglo+axis_offset+tic_len
+         WRITE (997,FMT='("closepath stroke")')
+         ! Tic label
+         WRITE (997,FMT='(I3," u ",1x,I3," u moveto")') i,jpjglo+axis_offset+tic_len+tic_label_space
+         WRITE (intStr, FMT='(I4)') i
+         WRITE (997,FMT='("(",(A),") dup stringwidth pop 2 div neg 0 rmoveto show")') TRIM(ADJUSTL(intStr))
+      END DO
+ 
+      ! Now draw the outline of each individual domain, nprocp should have been
+      ! set at the very beginning of the recursive partioning process.
+      WRITE (997,FMT='("0.7 setlinewidth")') ! Back to default width
+      WRITE (997,FMT='("0 0 0 setrgbcolor")')! and colour
+      DO iproc=1,nprocp
+         WRITE (997,FMT='("newpath")')
+         WRITE (997,FMT='("% Cursor initialization")')
+         WRITE (997,FMT='(I3," u ",1x,I3," u moveto %BL Corner")') pielb(iproc),pjelb(iproc)
+         WRITE (997,FMT='("% Drawing the rectangle")')
+         WRITE (997,FMT='(I3," u ",1x,I3," u lineto %BR Corner")') pieub(iproc),pjelb(iproc)
+         WRITE (997,FMT='(I3," u ",1x,I3," u lineto %TR Corner")') pieub(iproc),pjeub(iproc)
+         WRITE (997,FMT='(I3," u ",1x,I3," u lineto %TL Corner")') pielb(iproc),pjeub(iproc)
+         WRITE (997,FMT='("closepath stroke")')
+         WRITE (997,FMT='(I3," u ",1x,I3," u moveto")') &
+                    INT(0.5*(pieub(iproc)+pielb(iproc))), &
+                    INT(0.5*(pjeub(iproc)+pjelb(iproc)-1))
+         ! Write processor label, left justified
+         WRITE (intStr, FMT='(I4)') iproc-1
+         ! Use postscipt operator 'stringwidth' to get the approximate width of the
+         ! string containing the PE id and then adjust position so it is centred
+         ! about point we've just moveto'd. This doesn't attempt to deal with
+         ! vertical centering.
+         WRITE (997,FMT='("(",(A),") dup stringwidth pop 2 div neg 0 rmoveto show")') TRIM(ADJUSTL(intStr))
+      END DO
+      WRITE (997,FMT='("showpage")')
+      WRITE (997,FMT='("%%EOF")')
+      CLOSE(997)
+
+    END SUBROUTINE write_partition
+
+END MODULE partition_mod
