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
