Index: NEMO/branches/UKMO/dev_fix_mpp_prep/MPP_PREP/src/mpp_optimiz_zoom_nc.f90
===================================================================
--- NEMO/branches/UKMO/dev_fix_mpp_prep/MPP_PREP/src/mpp_optimiz_zoom_nc.f90	(revision 10264)
+++ NEMO/branches/UKMO/dev_fix_mpp_prep/MPP_PREP/src/mpp_optimiz_zoom_nc.f90	(revision 10265)
@@ -1,375 +1,535 @@
-PROGRAM mpp_optimiz_nc
- !!---------------------------------------------------------------------
- !!
- !!                       PROGRAM MPP_OPTIMIZ_NC
- !!                     ***********************
- !!
- !!  PURPOSE :
- !!  ---------
- !!              This program is build to optimize the domain beakdown into
- !!              subdomain for mpp computing.
- !!              Once the grid size, and the land/sea mask is known, it looks
- !!              for all the possibilities within a range of setting parameters
- !!              and determine the optimal.
- !!
- !!              Optimization is done with respect to the maximum number of
- !!              sea processors and to the maximum numbers of procs (jprocx)
- !!                     
- !!              Optional optimization can be performed takink into account
- !!              the maximum available processor memory ppmcal. This is
- !!              activated if jpmen =1
- !!
- !! history:
- !! --------
- !!       original  : 95-12 (Imbard M) for OPA8.1, CLIPPER
- !!       f90       : 03-06 (Molines JM), namelist as input
- !!                 : 05-05 (Molines JM), bathy in ncdf
- !!----------------------------------------------------------------------
- !! * modules used
-  USE netcdf
-
-  IMPLICIT NONE
-
-  INTEGER ::  jprocx=250   !: maximum number of proc. (Read from namelist)
-  INTEGER ::  jpmem=0      !: memory constraint (1) or no constraint (0)
-     !                     !  (use 1 with caution as the memory size of 
-     !                     !   the code lays on OPA 8.1 estimates ...)
-     !
-  INTEGER          ::  &
-       jpk    = 46  ,    & !: vertical levels (namelist)
-       jpiglo = 1442,    & !: I-size of the model (namelist)
-       jpjglo = 1021,    & !: J-size of the model (namelist)
-       jpidta = 1442,    & !: I-size of the data file (namelist)
-       jpjdta = 1021,   &  !: J-size of the data files (namelist)
-       nizoom = 1 ,     &  !: I zoom indicator (namelist)
-       njzoom = 1 ,     &  !: J zoom indicatori (namelist)
-       numnam = 4          !: logical unit for the namelist
-  NAMELIST /namspace/ jpk,jpiglo,jpjglo,jpidta,jpjdta,nizoom,njzoom
-  NAMELIST /namproc/ jprocx, jpmem
-
-  INTEGER ::  jpnix ,jpnjx  
+PROGRAM mpp_optimize
+   !!======================================================================
+   !!                     ***  PROGRAM  mpp_optimize  ***
+   !!=====================================================================
+   !!  ** Purpose : Propose possible domain decompositions for a given 
+   !!               bathymetric file, which is particularly intersting when
+   !!               we want to eliminate land-only domain. 
+   !!               All solution are proposed and written to output file.
+   !!               The ratio between the effective number of computed 
+   !!               point and the total number of points in the domain is 
+   !!               given and is probably a major criteria for choosing a 
+   !!               domain decomposition.
+   !!
+   !!  ** Method  : Use mpp_init like code for seting up the decomposition
+   !!               and evaluate the efficiency of the decomposition.
+   !! History
+   !!       original  : 95-12 (Imbard M) for OPA8.1, CLIPPER
+   !!       f90       : 03-06 (Molines JM), namelist as input
+   !!                 : 05-05 (Molines JM), bathy in ncdf
+   !!                 : 13-03 (Molines JM), Nemo-like coding and license.
+   !!                 : 18-10 (Mathiot  P), upgrade the NEMO 4.0
+   !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
+   !!   routines      : description
+   !!----------------------------------------------------------------------
+
+
+   !!----------------------------------------------------------------------
+   !! MPP-PREP, MEOM 2013
+   !! $Id$
+   !! Copyright (c) 2013, J.-M. Molines
+   !! Software governed by the CeCILL licence (Licence/MPP-PREPCeCILL.txt)
+   !!----------------------------------------------------------------------
+   USE netcdf
+
+   IMPLICIT NONE
+
+   INTEGER, PARAMETER :: jpreci=1 ,jprecj=1   !: overlap between processors
+
+   ! Namelist declaration and definition
+   ! -----------------------------------
+   INTEGER ::  nn_procmax  =250    !: maximum number of proc. (Read from namelist)
+   INTEGER ::  nn_procmin  = 1     !: maximum number of proc. (Read from namelist)
+   LOGICAL ::  ln_memchk = .FALSE. ! add a memory constraint if true (obsolete)
+   NAMELIST /namproc/ nn_procmax, nn_procmin, ln_memchk
+   !
+   INTEGER ::  nn_jpk = 46   !: vertical levels 
+   INTEGER ::  nn_izoom = 1  !: I zoom indicator
+   INTEGER ::  nn_jzoom = 1  !: J zoom indicator
+   NAMELIST /namspace/ nn_jpk, nn_izoom, nn_jzoom
+   !
+   ! Following variables are used only if ln_memchk=.true.
+   REAL(KIND=4) ::  required_memory, rppmpt !: not in namelist working array
+   REAL(KIND=4) ::  rn_ppmcal = 225000000. !: maximum memory of one processor for a 
+   !: given machine (in 8 byte words)
+   REAL(KIND=4) ::  rn_ppmin  = 0.4        !: minimum ratio to fill the memory
+   REAL(KIND=4) ::  rn_ppmax = 0.9         !: maximum ratio to fill the memory
+   NAMELIST /namparam/ rn_ppmcal, rn_ppmin, rn_ppmax
+   !
+   CHARACTER(LEN=80) :: cn_var='none'   !: Variable name of the bathymetry
+   CHARACTER(LEN=80) :: cn_x='x'        !: X dimension name
+   CHARACTER(LEN=80) :: cn_y='y'        !: Y dimension name
+   CHARACTER(LEN=80) :: cn_fbathy       !: File name of the netcdf bathymetry (namelist)
+   LOGICAL           :: ln_zps=.FALSE.  !: Logical flag for partial cells.
+   NAMELIST /namfile/ cn_fbathy, cn_var, cn_x, cn_y,  ln_zps
+   !
+   CHARACTER(LEN=80) :: cn_fovdta     !: root file name for keep output
+   NAMELIST /namkeep/ cn_fovdta
+   !
+   INTEGER            :: numnam = 4       ! logical unit for namelist
+   INTEGER            :: numout = 10        ! logical unit for output
+   INTEGER            :: npiglo, npjglo   ! domain size
+   INTEGER            :: npidta, npjdta   ! domain size
+
+   INTEGER            :: ji, jj, jni, jnj ! dummy loop index
+   INTEGER            :: ii, ij, jjc  ! dummy loop index
+   INTEGER            :: narg, iargc, ijarg      ! browsing command line
+
+   ! Decomposition related arrays (using same meaning than in NEMO)
+   INTEGER, DIMENSION(:,:), ALLOCATABLE :: ilci, ilcj ,iimppt, ijmppt
+   INTEGER, DIMENSION(:)  , ALLOCATABLE :: nlei_ocea, nldi_ocea
+   INTEGER, DIMENSION(:)  , ALLOCATABLE :: nlej_ocea, nldj_ocea
+   INTEGER, DIMENSION(:)  , ALLOCATABLE :: nlei_land, nldi_land
+   INTEGER, DIMENSION(:)  , ALLOCATABLE :: nlej_land, nldj_land
+   INTEGER                              :: nimpp, njmpp
+   INTEGER                              :: nreci, nrecj
+   INTEGER                              :: ili, ilj
+   INTEGER                              :: jarea, iarea, iarea0
+   INTEGER                              :: iresti, irestj
+   !
+   INTEGER :: ioce, isurf             !: number of ocean points cumulated, per_proc
+   INTEGER :: ioce_opt                !: number of ocean points cumulated for optimal case
+   INTEGER :: nland, nocea, nvalid    !: number of land, ocean, memory_valid  procs 
+   INTEGER :: nland_opt               !: optimal number of land procs
+   INTEGER :: ii1, ii2, ij1, ij2      !: limit of subdomain in global domain
+   INTEGER :: jpimax,     jpjmax            !: size of sub domain
+   INTEGER :: jpimax_opt, jpjmax_opt        !: size of sub domain for optimal case
+   INTEGER :: inf10,     inf30,     inf50      !: 
+   INTEGER :: inf10_opt, inf30_opt, inf50_opt  !:  in optimal case
+   INTEGER :: npni_opt, npnj_opt      !: optimal domain decomposition
+
+   INTEGER :: iminproci, imaxproci    !: limits of the processor loop
+   INTEGER :: iminprocj, imaxprocj    !: can be reduded to  nkeepi, nkeepj
+
+   ! Saving criteria
+   REAL(KIND=4) :: ratio_min=99999.   !: keep only decomposition with ration less than ratio_min
+   INTEGER      :: nocea_min = 1      !: minimum number of ocean procs for saving
+   INTEGER      :: nmodulo = 1        !: Only keep solution multiple of nmodulo
+   LOGICAL      :: ll_criteria=.TRUE. !:
+   !
+   REAL(KIND=4)                              ::  oce_cover
+   REAL(KIND=4)                              ::  oce_cover_min,     oce_cover_max,     ratio
+   REAL(KIND=4)                              ::  oce_cover_min_opt, oce_cover_max_opt, ratio_opt
+   REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE ::  tmask     ! npiglo x npjglo
+   REAL(KIND=4), DIMENSION(:,:), ALLOCATABLE ::  bathy     ! npidta x npjdta
+
+   ! CDF stuff
+   INTEGER :: ncid, istatus, id
+   LOGICAL ::  ll_good = .FALSE.
+
+   CHARACTER(LEN=80) :: cf_namlist='namelist'
+   CHARACTER(LEN=80) :: cf_out='processor.layout'
+   CHARACTER(LEN=80) :: cdum                       ! dummy character variable
+
+   ! Keep stuff      
+   LOGICAL ::  ll_keep = .FALSE.
+   INTEGER :: nkeepi, nkeepj          !: for option -keep : the retained decomposition
   !
-  INTEGER,PARAMETER :: jpreci=1 ,jprecj=1
-  !
-  ! Following variables are used only if jpmem=1
-  REAL(KIND=4) ::  ppmpt ,   &
-       ppmcal = 225000000., &  !: maximum memory of one processor for a given machine (in 8 byte words)
-       ppmin  = 0.4,         & !: minimum ratio to fill the memory
-       ppmax  = 0.9            !: maximum ration to fill the memory
-  ! Aleph
-  !     PARAMETER(ppmcal= 16000000.)
-  !Brodie
-  !     PARAMETER(ppmcal=250000000.)
-  ! Uqbar
-  !     PARAMETER(ppmcal=3750000000.)
-  ! Zahir
-  !     PARAMETER(ppmcal=225000000.)
-
-  CHARACTER(LEN=80) :: cbathy, &       !: File name of the netcdf bathymetry (namelist)
-      &                clvar           !: Variable name in netcdf for the bathy to be read
-  LOGICAL ::  ln_zps=.false.           !: Logical flag for partial cells.
-  NAMELIST /namfile/ cbathy, ln_zps
-  NAMELIST /namparam/ ppmcal, ppmin, ppmax
-  !
-  INTEGER :: iumout = 1
-  INTEGER :: ji,jj,jn,jni,jnj,jni2,jnj2
-  INTEGER :: iumbat,ifreq,il1,il2
-  INTEGER :: ii,iim,ij,ijm,imoy,iost,iresti,irestj,isurf,ivide
-  INTEGER :: iilb,ijlb,ireci,irecj,in
-  INTEGER :: ipi,ipj
-  INTEGER :: inf10,inf30,inf50,iptx,isw
-  INTEGER :: iii,iij,iiii,iijj,iimoy,iinf10,iinf30,iinf50
-  !
-  INTEGER,DIMENSION(:,:),ALLOCATABLE     ::  ibathy    ! jpidta -jpjdta
-  INTEGER,DIMENSION(:,:),ALLOCATABLE     ::  ippdi, ippdj ,iidom, ijdom
-  !
-  REAL(KIND=4)                           ::  zmin,zmax,zper,zmem
-  REAL(KIND=4)                           ::  zzmin,zzmax,zperx
-  REAL(KIND=4),DIMENSION(:,:),ALLOCATABLE  ::  zmask ,&  ! jpiglo -jpjglo
-      &                                        zdta      ! jpidta -jpjdta
-
- ! CDF stuff
-  INTEGER :: ncid, ivarid, istatus
-  LOGICAL ::  llbon=.false.
-  !
-  ! 0. Initialisation
-  ! -----------------
-  OPEN(numnam,FILE='namelist')
-  REWIND(numnam)
-  READ(numnam,namspace)
-
-  REWIND(numnam)
-  READ(numnam,namfile)
-
-  REWIND(numnam)
-  READ(numnam,namparam)
-
-  REWIND(numnam)
-  READ(numnam,namproc)
-
-  ! estimated  code size expressed in number of 3D arrays (valid for OPA8.1)
-  ppmpt = 55.+73./jpk
-  jpnix = jprocx ; jpnjx=jprocx
-
-  ALLOCATE ( ibathy(jpidta,jpjdta), zmask(jpiglo,jpjglo),zdta(jpidta,jpjdta) )
-  ALLOCATE (ippdi(jpnix,jpnjx), ippdj(jpnix,jpnjx) )
-  ALLOCATE (iidom(jpnix,jpnjx), ijdom(jpnix,jpnjx) )
-
-  OPEN(iumout,FILE='processor.layout')
-  WRITE(iumout,*)
-  WRITE(iumout,*) ' optimisation de la partition'
-  WRITE(iumout,*) ' ----------------------------'
-  WRITE(iumout,*)
-  !
-  ! * Read cdf bathy file
-  !
-         IF ( ln_zps ) THEN        ! partial steps
-            clvar = 'Bathymetry'
-         ELSE 
-            clvar = 'Bathy_level'  ! full steps
+
+   !!----------------------------------------------------------------------
+   narg=iargc()
+   ijarg=1
+   IF ( narg == 0 ) THEN
+      PRINT *,' try mpp_optimize -h for instructions !'
+      STOP
+   ENDIF
+   !
+   DO WHILE ( ijarg <= narg )
+      CALL getarg(ijarg,cdum) ; ijarg=ijarg+1
+      SELECT CASE ( cdum )
+      CASE ('-h') 
+         PRINT *,'  usage : mpp_optimize [ -h ]  [-keep jpni jpnj] [ -o file out ] '
+         PRINT *,'               [ -modulo val ] [-r ratio] [-minocean procs] -n namelist'
+         PRINT *,'      '
+         PRINT *,'     PURPOSE :'
+         PRINT *,'         This program is build to optimize the domain beakdown into'
+         PRINT *,'         subdomain for mpp computing.'
+         PRINT *,'         Once the grid size, and the land/sea mask is known, it looks'
+         PRINT *,'         for all the possibilities within a range of setting parameters'
+         PRINT *,'         and determine the optimal.'
+         PRINT *,''
+         PRINT *,'         Optimization is done with respect to the maximum number of'
+         PRINT *,'         sea processors and to the maximum numbers of procs (nn_procmax)'
+         PRINT *,'                '
+         PRINT *,'         Optional optimization can be performed taking into account'
+         PRINT *,'         the maximum available processor memory rn_ppmcal. This is'
+         PRINT *,'         activated if ln_memchk is set true in the namelist'
+         PRINT *,'      '
+         PRINT *,'         Additional criteria can be given on the command line to reduce'
+         PRINT *,'         the amount of possible choices.'
+         PRINT *,'      '
+         PRINT *,'     ARGUMENTS :'
+         PRINT *,'         -n namelist : indicate the name of the namelist to use'
+         PRINT *,'      '
+         PRINT *,'     OPTIONS :'
+         PRINT *,'         -h : print this help message'
+         PRINT *,'         -keep jpni jpnj : print a file suitable for plotting,'
+         PRINT *,'                 corresponding to the given decomposition'
+         PRINT *,'         -o output file : give the name of the output file'
+         PRINT *,'                 default is ',TRIM(cf_out)
+         PRINT *,'         -modulo val : only retain decomposition whose total number'
+         PRINT *,'                 of util processors (sea) are a multiple of val'
+         PRINT *,'         -r ratio : only retain decomposition with a ratio computed/global'
+         PRINT *,'                 less or equal to the given ratio'
+         PRINT *,'         -minocean procs : only retain decomposition with a number of '
+         PRINT *,'                 ocean procs greater of equal to procs'
+         PRINT *,'      '
+         PRINT *,'     REQUIRED FILES :'
+         PRINT *,'       A bathymetric file and an ad-hoc namelist are required.'
+         PRINT *,'       The file name of the bathymetry is specified in the namelist'
+         PRINT *,'      '
+         PRINT *,'     OUTPUT : '
+         PRINT *,'       ',TRIM(cf_out),' : an ascii file with all found possibilities'
+         PRINT *,'      '
+         STOP
+      CASE ('-n' )
+         CALL getarg(ijarg,cf_namlist) ; ijarg=ijarg+1
+      CASE ('-o' )
+         CALL getarg(ijarg,cf_out) ; ijarg=ijarg+1
+      CASE ('-keep' )
+         ll_keep=.TRUE.
+         CALL getarg(ijarg,cdum) ; ijarg=ijarg+1 ; READ( cdum,*) nkeepi
+         CALL getarg(ijarg,cdum) ; ijarg=ijarg+1 ; READ( cdum,*) nkeepj
+      CASE ('-modulo' )
+         CALL getarg(ijarg,cdum) ; ijarg=ijarg+1 ; READ( cdum,*) nmodulo
+      CASE ('-r' )
+         CALL getarg(ijarg,cdum) ; ijarg=ijarg+1 ; READ( cdum,*) ratio_min
+      CASE ('-minocean' )
+         CALL getarg(ijarg,cdum) ; ijarg=ijarg+1 ; READ( cdum,*) nocea_min
+      END SELECT
+   ENDDO
+
+   ! Open and read the namelist
+   OPEN(numnam,FILE=cf_namlist)
+   REWIND(numnam)
+   READ(numnam,namspace)
+
+   REWIND(numnam)
+   READ(numnam,namfile)
+
+   REWIND(numnam)
+   READ(numnam,namparam)
+
+   REWIND(numnam)
+   READ(numnam,namproc)
+
+   REWIND(numnam)
+   READ(numnam,namkeep)  ! only used for -keep option but still ...
+   CLOSE(numnam)
+
+   ! estimated code size expressed in number of 3D arrays (valid for OPA8.1) to be tuned for OPA9.0/Nemo
+   rppmpt = 55.+73./nn_jpk
+
+   ! Open bathy file an allocate required memory
+   INQUIRE( FILE=cn_fbathy, EXIST=ll_good )
+   IF( ll_good ) THEN
+      istatus = NF90_OPEN(cn_fbathy, NF90_NOWRITE, ncid)
+      istatus = NF90_INQ_DIMID(ncid, cn_x, id) ; istatus = NF90_INQUIRE_DIMENSION(ncid, id, len=npiglo)
+      istatus = NF90_INQ_DIMID(ncid, cn_y, id) ; istatus = NF90_INQUIRE_DIMENSION(ncid, id, len=npjglo)
+      npidta  = npiglo ; npjdta=npjglo
+   ELSE
+      PRINT *,' File missing : ', TRIM(cn_fbathy)
+      STOP 42
+   ENDIF
+
+   ALLOCATE (tmask(npiglo,npjglo), bathy(npidta,npjdta) )
+   ALLOCATE (ilci(nn_procmax,nn_procmax), ilcj(nn_procmax,nn_procmax) )
+   ALLOCATE (iimppt(nn_procmax,nn_procmax), ijmppt(nn_procmax,nn_procmax) )
+
+   ! Open output file for results
+   IF ( ll_keep ) THEN
+      nn_procmax = nkeepi*nkeepj  ! reduce nn_procmax
+      ! File will be open later
+   ELSE
+      OPEN(numout,FILE=cf_out)
+      WRITE(numout,*)
+      WRITE(numout,*) ' Domain decomposition optimization '
+      WRITE(numout,*) ' ----------------------------------'
+      WRITE(numout,*)
+   ENDIF
+   !
+   ALLOCATE ( nlei_ocea(nn_procmax), nldi_ocea(nn_procmax), nlej_ocea(nn_procmax), nldj_ocea(nn_procmax) )
+   ALLOCATE ( nlei_land(nn_procmax), nldi_land(nn_procmax), nlej_land(nn_procmax), nldj_land(nn_procmax) )
+   !
+   ! Read cdf bathy file
+   IF ( cn_var == 'none' ) THEN  ! automatic detection of variable name according to partial step
+      IF ( ln_zps ) THEN           ! partial steps
+         cn_var = 'Bathymetry'
+      ELSE 
+         cn_var = 'Bathy_level'    ! full steps
+      ENDIF
+   ENDIF
+   PRINT *,''
+   PRINT *,' ocean/land file used is: ', TRIM(cn_fbathy)
+   PRINT *,' variable used to find ocean domain is: ', TRIM(cn_var)
+   PRINT *,' Dimensions (jpi x jpj) are: ',npiglo,'x',npjglo
+   PRINT *,''
+
+   istatus = NF90_INQ_VARID (ncid, cn_var, id)
+   istatus = NF90_GET_VAR   (ncid, id,   bathy)
+   istatus = NF90_CLOSE     (ncid)
+   !
+   ! Building the mask ( eventually on a smaller domain than the bathy)
+   tmask(:,:) = bathy(nn_izoom:nn_izoom+npiglo -1,  nn_jzoom:nn_jzoom+npjglo -1)
+
+   WHERE ( tmask > 0 ) 
+      tmask = 1.
+   ELSEWHERE
+      tmask = 0.
+   ENDWHERE
+
+   !  Main loop on processors
+   ! ------------------------
+   ! initialization of working variables
+   npni_opt=1       ; npnj_opt=1
+   jpimax_opt=npiglo ; jpjmax_opt=npjglo
+   nland_opt=0   
+   ioce_opt=0
+   oce_cover_min_opt=0. ; oce_cover_max_opt=0.
+   inf10_opt=0 ; inf30_opt=0 ; inf50_opt=0
+   ratio_opt=1.
+
+   nvalid=0       ! counter for valid case ( ln_memchk true )
+   IF ( ll_keep ) THEN
+      iminproci = nkeepi    ; imaxproci = iminproci
+      iminprocj = nkeepj    ; imaxprocj = iminprocj
+   ELSE
+      iminproci = 1    ; imaxproci = MIN( nn_procmax, npiglo )
+      iminprocj = 1    ; imaxprocj = MIN( nn_procmax, npjglo )
+   ENDIF
+
+   ! loop on all decomposition a priori
+   PRINT *, 'Loop over all the decomposition (can take a while) ...'
+   PRINT *, ''
+   DO jni=iminproci, imaxproci
+      DO jnj=iminprocj, imaxprocj
+         ! Limitation of the maxumun number of PE's
+         IF ( jni*jnj <=  nn_procmax .AND. jni*jnj >= nn_procmin )  THEN
+            !
+            !  1. Dimension arrays for subdomains
+            ! -----------------------------------
+            !
+            ! Partition : size of sub-domain 
+            jpimax=(npiglo-2*jpreci + (jni-1))/jni + 2*jpreci
+            jpjmax=(npjglo-2*jprecj + (jnj-1))/jnj + 2*jprecj
+            !
+            ! Memory optimization ?
+            IF ( ln_memchk ) THEN
+               required_memory=rppmpt*jpimax*jpjmax*nn_jpk
+               IF( required_memory > rn_ppmcal ) EXIT
+               IF( required_memory > rn_ppmax*rn_ppmcal .OR. required_memory < rn_ppmin*rn_ppmcal) EXIT
+            ENDIF
+            nvalid=nvalid+1
+            !
+            ! Position of each sub domain   (jni x jni in total )
+            nreci  = 2*jpreci                      ; nrecj  = 2*jprecj
+            iresti = 1 + MOD ( npiglo - nreci - 1 , jni )  ; irestj = 1 + MOD ( npjglo - nrecj - 1 , jnj )
+            !
+            ! 
+            ilci(       1:iresti, 1:jnj) = jpimax
+            ilci(iresti+1:jni   , 1:jnj) = jpimax-1
+
+            ilcj(1:jni,       1:irestj) = jpjmax
+            ilcj(1:jni,irestj+1:jnj   ) = jpjmax-1
+
+            !  2. Index arrays for subdomains
+            ! -------------------------------
+            iimppt(1:jni, 1:jnj) =  1
+            ijmppt(1:jni, 1:jnj) =  1
+            IF( jni > 1 ) THEN
+               DO jj=1,jnj
+                  DO ji=2,jni
+                     iimppt(ji,jj)= iimppt(ji-1,jj) + ilci(ji-1,jj) - nreci
+                  END DO
+               END DO
+            ENDIF
+
+            IF( jnj > 1 ) THEN
+               DO jj=2,jnj
+                  DO ji=1,jni
+                     ijmppt(ji,jj)= ijmppt(ji,jj-1) + ilcj(ji,jj-1) - nrecj
+                  END DO
+               END DO
+            ENDIF
+            !
+            ! Loop on each subdomain to look for land proportion
+            nland = 0
+            nocea = 0
+            ioce  = 0
+            oce_cover_min = 1.e+20
+            oce_cover_max = -1.e+20
+            inf10=0
+            inf30=0
+            inf50=0
+            !
+            ! 3. Subdomain description in the Regular Case
+            ! --------------------------------------------
+            !
+            DO jarea = 1, jni*jnj
+                  iarea0 = jarea - 1
+                  ii = 1 + MOD(iarea0,jni)
+                  ij = 1 +     iarea0/jni
+                  ili = ilci(ii,ij)
+                  ilj = ilcj(ii,ij)
+
+                  isurf = 0
+                  ! loop on inner point of sub-domain
+                  DO jj=1, ilj
+                     DO  ji=1, ili
+                        IF( tmask(ji + iimppt(ii,ij) - 1, jj + ijmppt(ii,ij) - 1) == 1 ) isurf=isurf+1
+                     END DO
+                  END DO
+
+                  nimpp = iimppt(ii,ij)
+                  njmpp = ijmppt(ii,ij)
+                  ii1   = nimpp+jpreci      ; ii2 = nimpp+ili-1 -jpreci
+                  ij1   = njmpp+jprecj      ; ij2 = njmpp+ilj-1 -jprecj
+                  IF ( isurf == 0 ) THEN
+                     nland = nland+1
+                     nldi_land(nland) = ii1
+                     nlei_land(nland) = ii2
+                     nldj_land(nland) = ij1
+                     nlej_land(nland) = ij2
+                  ELSE
+                     nocea = nocea+1
+                     ioce  = ioce + isurf
+                     nldi_ocea(nocea) = ii1
+                     nlei_ocea(nocea) = ii2
+                     nldj_ocea(nocea) = ij1
+                     nlej_ocea(nocea) = ij2
+                  ENDIF
+
+                  ! ratio of wet points over total number of point per proc.
+                  oce_cover = float(isurf)/float(jpimax*jpjmax)
+
+                  IF(oce_cover_min > oce_cover .AND. isurf /= 0) oce_cover_min=oce_cover
+                  IF(oce_cover_max < oce_cover .AND. isurf /= 0) oce_cover_max=oce_cover
+                  IF(oce_cover     < 0.1       .AND. isurf /= 0) inf10=inf10+1
+                  IF(oce_cover     < 0.3       .AND. isurf /= 0) inf30=inf30+1
+                  IF(oce_cover     < 0.5       .AND. isurf /= 0) inf50=inf50+1
+                  !
+               !END DO  ! loop on processors
+            END DO     ! loop on processors
+            ! 
+            ratio=float(nocea)*float(jpimax*jpjmax)/float(npiglo*npjglo)
+
+            ! criteria for printing results
+            ll_criteria = ( ( MOD ( nocea, nmodulo ) == 0 ) .AND. &
+                 &          ( ratio <= ratio_min          ) .AND. &
+                 &          ( nocea >= nocea_min           )  )
+            IF ( ll_keep ) THEN   ! the loop in done only once !
+               WRITE(cdum,'(a,"-",i3.3,"x",i3.3,"_",i4.4)') TRIM(cn_fovdta), nkeepi, nkeepj, nocea
+               OPEN(numout, file=cdum )
+               WRITE(numout,'("# ocean ",i5)') nocea
+               DO jjc=1, nocea
+                  WRITE(numout,'("#",i5)') jjc
+                  WRITE(numout,'(2i5)') nldi_ocea(jjc)-1+nn_izoom-1, nldj_ocea(jjc)-1+nn_jzoom -1
+                  WRITE(numout,'(2i5)') nlei_ocea(jjc)+1+nn_izoom-1, nldj_ocea(jjc)-1+nn_jzoom -1
+                  WRITE(numout,'(2i5)') nlei_ocea(jjc)+1+nn_izoom-1, nlej_ocea(jjc)+1+nn_jzoom -1
+                  WRITE(numout,'(2i5)') nldi_ocea(jjc)-1+nn_izoom-1, nlej_ocea(jjc)+1+nn_jzoom -1
+                  WRITE(numout,'(2i5)') nldi_ocea(jjc)-1+nn_izoom-1, nldj_ocea(jjc)-1+nn_jzoom -1
+                  WRITE(numout,'(2i5)') 9999, 9999
+               ENDDO
+               !
+               WRITE(numout,'("# land ",i5)') nland
+               DO jjc=1, nland
+                  WRITE(numout,'("# land ",i5)') jjc
+                  WRITE(numout,'(2i5)') nldi_land(jjc)-1+nn_izoom-1, nldj_land(jjc)-1+nn_jzoom -1
+                  WRITE(numout,'(2i5)') nlei_land(jjc)+1+nn_izoom-1, nldj_land(jjc)-1+nn_jzoom -1
+                  WRITE(numout,'(2i5)') nlei_land(jjc)+1+nn_izoom-1, nlej_land(jjc)+1+nn_jzoom -1
+                  WRITE(numout,'(2i5)') nldi_land(jjc)-1+nn_izoom-1, nlej_land(jjc)+1+nn_jzoom -1
+                  WRITE(numout,'(2i5)') nldi_land(jjc)-1+nn_izoom-1, nldj_land(jjc)-1+nn_jzoom -1
+                  WRITE(numout,'(2i5)') nlei_land(jjc)+1+nn_izoom-1, nlej_land(jjc)+1+nn_jzoom -1
+                  WRITE(numout,'(2i5)') nldi_land(jjc)-1+nn_izoom-1, nlej_land(jjc)+1+nn_jzoom -1
+                  WRITE(numout,'(2i5)') nlei_land(jjc)+1+nn_izoom-1, nldj_land(jjc)-1+nn_jzoom -1
+                  WRITE(numout,'(2i5)') 9999, 9999
+               ENDDO
+               !
+            ELSE
+               IF ( ll_criteria ) THEN
+                  WRITE(numout,*) ' iresti=',iresti,' irestj=',irestj
+                  WRITE(numout,*) '--> Total number of domains ',jni*jnj
+                  WRITE(numout,*) ' '
+                  WRITE(numout,*) ' jpni=',jni ,' jpnj=',jnj
+                  WRITE(numout,*) ' jpi= ',jpimax ,' jpj= ',jpjmax
+                  WRITE(numout,*) ' Number of ocean processors       ', nocea
+                  WRITE(numout,*) ' Number of land processors        ', nland
+                  WRITE(numout,*) ' Mean ocean coverage per domain   ', float(ioce)/float(nocea)/float(jpimax*jpjmax)
+                  WRITE(numout,*) ' Minimum ocean coverage           ', oce_cover_min
+                  WRITE(numout,*) ' Maximum ocean coverage           ', oce_cover_max
+                  WRITE(numout,*) ' nb of proc with coverage         < 10 % ', inf10
+                  WRITE(numout,*) ' nb of proc with coverage 10 < nb < 30 % ', inf30 - inf10
+                  WRITE(numout,*) ' nb of proc with coverage 30 < nb < 50 % ', inf50 - inf30
+                  WRITE(numout,*) ' Number of computed points        ', nocea*jpimax*jpjmax
+                  WRITE(numout,*) ' Overhead of computed points      ', nocea*jpimax*jpjmax-npiglo*npjglo
+                  WRITE(numout,*) ' % sup (computed / global)        ', ratio
+                  WRITE(numout,*)
+               ENDIF   ! note that indication of optimum does not take modulo into account (for information)
+               ! 
+               ! Look for optimum 
+               IF( nland > nland_opt ) THEN
+                  npni_opt          = jni
+                  npnj_opt          = jnj
+                  jpimax_opt           = jpimax
+                  jpjmax_opt           = jpjmax
+                  nland_opt         = nland
+                  ioce_opt          = ioce
+                  oce_cover_min_opt = oce_cover_min
+                  oce_cover_max_opt = oce_cover_max
+                  inf10_opt         = inf10
+                  inf30_opt         = inf30
+                  inf50_opt         = inf50
+                  ratio_opt         = ratio
+               ELSE IF( nland == nland_opt .AND. ratio_opt < ratio) THEN
+                  npni_opt          = jni
+                  npnj_opt          = jnj
+                  jpimax_opt           = jpimax
+                  jpjmax_opt           = jpjmax
+                  ioce_opt          = ioce
+                  oce_cover_min_opt = oce_cover_min
+                  oce_cover_max_opt = oce_cover_max
+                  inf10_opt         = inf10
+                  inf30_opt         = inf30
+                  inf50_opt         = inf50
+                  ratio_opt         = ratio
+               ENDIF
+            ENDIF
          ENDIF
-
-         INQUIRE( FILE=cbathy, EXIST=llbon )
-      IF( llbon ) THEN
-            istatus=NF90_OPEN(cbathy,NF90_NOWRITE,ncid)
-            istatus=NF90_INQ_VARID(ncid,clvar,ivarid)
-            istatus=NF90_GET_VAR(ncid,ivarid,zdta)
-            istatus=NF90_CLOSE(ncid)
-      ELSE
-          PRINT *,' File missing : ', trim(cbathy)
-          STOP
+      END DO
+   END DO
+   !
+   ! print optimal result
+   IF ( .NOT. ll_keep ) THEN
+      IF ( nvalid == 0 ) THEN
+         WRITE(numout,*) ' no possible choice ...'
+         WRITE(numout,*)
+         WRITE(numout,*) 'insufficient number of processors for the available memory'
+         STOP 
       ENDIF
-  ibathy(:,:)=zdta(:,:)
-
-  !
-  ! Building the mask
-  DO jj=1,jpjglo
-     DO ji=1,jpiglo
-        zmask(ji,jj) = float(ibathy(ji+nizoom - 1,jj+njzoom -1))
-     END DO
-  END DO
-
-  DO jj=1,jpjglo
-     DO ji=1,jpiglo
-        zmask(ji,jj)=  min(REAL(1.,kind=4),max(REAL(0.,kind=4),zmask(ji,jj)))  ! Old vector coding rule ...
-     END DO
-  END DO
-  !
-  !  Main loop on processors
-  ! ------------------------
-  iii=1 ; iij=1
-  iiii=jpiglo ; iijj=jpjglo
-  iptx=0
-  iimoy=0
-  zzmin=0. ; zzmax=0.
-  iinf10=0 ; iinf30=0 ; iinf50=0
-  zperx=1.
-  in=0
-  DO jni=1,jpnix
-     DO jnj=1,jpnjx
-        !
-        ! Limitation ob the maxumun number of PE's
-        IF(jni*jnj >  jprocx) goto 1000
-        !
-        ! Partition
-        ipi=(jpiglo-2*jpreci + (jni-1))/jni + 2*jpreci
-        ipj=(jpjglo-2*jprecj + (jnj-1))/jnj + 2*jprecj
-        !
-        ! Memory optimization ?
-        isw=0
-        zmem=ppmpt*ipi*ipj*jpk
-        IF(zmem > ppmcal) go to 1000
-        IF(jpmem == 1) THEN
-           IF(zmem.GT.ppmax*ppmcal.OR.zmem.LT.ppmin*ppmcal) isw=1
-        ENDIF
-        IF(isw.EQ.1) go to 1000
-        in=in+1
-        !
-        WRITE(iumout,*) '--> nombre de processeurs ',jni*jnj
-        WRITE(iumout,*) ' '
-        WRITE(iumout,*) " jpni=",jni ," jpnj=",jnj
-        WRITE(iumout,*) " jpi= ",ipi ," jpj= ",ipj
-        zper=(jni*jnj*ipi*ipj)/float(jpiglo*jpjglo)
-        WRITE(iumout,*) " rapport jpnij*domain/global domain ",zper
-        !
-        ! Coin en bas a gauche de chaque processeur
-        !
-        iilb=1
-        ijlb=1
-        ireci=2*jpreci
-        irecj=2*jprecj
-        iresti = MOD ( jpiglo - ireci , jni )
-        irestj = MOD ( jpjglo - irecj , jnj )
-        !
-        IF (iresti.EQ.0) iresti = jni
-        DO jj=1,jnj
-           DO ji=1,iresti
-              ippdi(ji,jj) = ipi
-           END DO
-           DO ji=iresti+1,jni
-              ippdi(ji,jj) = ipi -1
-           END DO
-        END DO
-        IF (irestj.EQ.0) irestj = jnj
-        DO ji=1,jni
-           DO jj=1,irestj
-              ippdj(ji,jj) = ipj
-           END DO
-           DO jj=irestj+1,jnj
-              ippdj(ji,jj) = ipj -1
-           END DO
-        END DO
-        DO jj=1,jnj
-           DO ji=1,jni
-              iidom(ji,jj)=iilb
-              ijdom(ji,jj)=ijlb
-           END DO
-        END DO
-        WRITE(iumout,*) " iresti=",iresti," irestj=",irestj
-        !
-        !  2. Boucle sur les processeurs
-        ! ------------------------------
-        !
-        ivide=0
-        imoy=0
-        zmin=1.e+20
-        zmax=-1.e+20
-        inf10=0
-        inf30=0
-        inf50=0
-        !
-        DO jni2=1,jni
-           DO jnj2=1,jnj
-
-              IF(jni.GT.1)THEN
-                 DO jj=1,jnj
-                    DO ji=2,jni
-                       iidom(ji,jj)=iidom(ji-1,jj)+ippdi(ji-1,jj)-ireci
-                    END DO
-                 END DO
-                 iilb=iidom(jni2,jnj2)
-              ENDIF
-              IF(jnj.GT.1)THEN
-                 DO jj=2,jnj
-                    DO ji=1,jni
-                       ijdom(ji,jj)=ijdom(ji,jj-1)+ippdj(ji,jj-1)-irecj
-                    END DO
-                 END DO
-                 ijlb=ijdom(jni2,jnj2)
-              ENDIF
-
-              ! Check wet points over the entire domain to preserve the MPI communication stencil
-              isurf=0
-              DO jj=1,ippdj(jni2,jnj2)
-                 DO  ji=1,ippdi(jni2,jnj2)
-                    IF(zmask(ji+iilb-1,jj+ijlb-1).EQ.1.) isurf=isurf+1
-                 END DO
-              END DO
-
-              IF(isurf.EQ.0) THEN
-                 ivide=ivide+1
-              ELSE
-                 imoy=imoy+isurf
-              ENDIF
-              zper=float(isurf)/float(ipi*ipj)
-              IF(zmin.GT.zper.AND.isurf.NE.0) zmin=zper
-              IF(zmax.LT.zper.AND.isurf.NE.0) zmax=zper
-              IF(zper.LT.0.1.AND.isurf.NE.0) inf10=inf10+1
-              IF(zper.LT.0.3.AND.isurf.NE.0) inf30=inf30+1
-              IF(zper.LT.0.5.AND.isurf.NE.0) inf50=inf50+1
-              !
-              !
-              ! 3. Fin de boucle sur les processeurs, impression
-              ! ------------------------------------------------
-              !
-           END DO
-        END DO
-        WRITE(iumout,*) ' nombre de processeurs       ',jni*jnj
-        WRITE(iumout,*) ' nombre de processeurs mer   ',jni*jnj-ivide
-        WRITE(iumout,*) ' nombre de processeurs terre ',ivide
-        WRITE(iumout,*) ' moyenne de recouvrement     ',float(imoy)/float(jni*jnj-ivide)/float(ipi*ipj)
-        WRITE(iumout,*) ' minimum de recouvrement     ',zmin
-        WRITE(iumout,*) ' maximum de recouvrement     ',zmax
-        WRITE(iumout,*) ' nb de p recouvrement < 10 % ',inf10
-        WRITE(iumout,*) ' nb de p      10 < nb < 30 % ',inf30-inf10
-        WRITE(iumout,*) ' nb de p      30 < nb < 50 % ',inf50-inf10 -inf30
-        WRITE(iumout,*) ' nombre de points integres   ', (jni*jnj-ivide)*ipi*ipj
-        WRITE(iumout,*) ' nbr de pts supplementaires  ', (jni*jnj-ivide)*ipi*ipj-jpiglo*jpjglo
-        zper=float((jni*jnj-ivide))*float(ipi*ipj)/float(jpiglo*jpjglo)
-        WRITE(iumout,*) ' % sup                       ',zper
-        WRITE(iumout,*)
-        ! 
-        ! 
-        ! 4. Recherche de l optimum
-        ! -------------------------
-        !
-        IF(ivide.GT.iptx) THEN
-           iii=jni
-           iij=jnj
-           iiii=ipi
-           iijj=ipj
-           iptx=ivide
-           iimoy=imoy
-           zzmin=zmin
-           zzmax=zmax
-           iinf10=inf10
-           iinf30=inf30
-           iinf50=inf50
-           zperx=zper
-        ELSE IF(ivide.EQ.iptx.AND.zperx.LT.zper) THEN
-           iii=jni
-           iij=jnj
-           iiii=ipi
-           iijj=ipj
-           iimoy=imoy
-           zzmin=zmin
-           zzmax=zmax
-           iinf10=inf10
-           iinf30=inf30
-           iinf50=inf50
-           zperx=zper
-        ENDIF
-        !
-        ! 5. Fin de boucle sur le nombre de processeurs
-        ! ---------------------------------------------
-        !
-      1000 continue
-     END DO
-  END DO
-  !
-  !
-  ! 6. Affichage resultat
-  ! ---------------------
-  !
-  IF(in.EQ.0) THEN
-     WRITE(iumout,*) ' le choix n'' a pas pu etre fait '
-     WRITE(iumout,*)
-     WRITE(iumout,*) 'le nombre de processeurs maximum est insuffisant'
-     STOP 
-  ENDIF
-  WRITE(iumout,*) ' choix optimum'
-  WRITE(iumout,*) ' ============='
-  WRITE(iumout,*) 
-  WRITE(iumout,*) '--> nombre de processeurs ',iii*iij
-  WRITE(iumout,*) ' '
-  WRITE(iumout,*) " jpni=",iii ," jpnj=",iij
-  WRITE(iumout,*) " jpi= ",iiii ," jpj= ",iijj
-  WRITE(iumout,*) 
-  WRITE(iumout,*) ' nombre de processeurs mer   ',iii*iij-iptx
-  WRITE(iumout,*) ' nombre de processeurs terre ',iptx
-  WRITE(iumout,*) ' moyenne de recouvrement     ',float(iimoy)/float(iii*iij-iptx)/float(iiii*iijj)
-  WRITE(iumout,*) ' minimum de recouvrement     ',zzmin
-  WRITE(iumout,*) ' maximum de recouvrement     ',zzmax
-  WRITE(iumout,*) ' nb de p recouvrement < 10 % ',iinf10
-  WRITE(iumout,*) ' nb de p      10 < nb < 30 % ',iinf30-iinf10
-  WRITE(iumout,*) ' nb de p      30 < nb < 50 % ',iinf50-iinf10 -iinf30
-  WRITE(iumout,*) ' nombre de points integres   ', (iii*iij-iptx)*iiii*iijj
-  WRITE(iumout,*) ' nbr de pts supplementaires  ', (iii*iij-iptx)*iiii*iijj-jpiglo*jpjglo
-  WRITE(iumout,*) ' % sup                       ',zperx
-  WRITE(iumout,*)
-  CLOSE(iumout)
-  !
-  !
-  !
-  STOP
-END PROGRAM mpp_optimiz_nc
+
+      WRITE(numout,*) ' Optimal choice'
+      WRITE(numout,*) ' =============='
+      WRITE(numout,*) 
+      WRITE(numout,*) '--> Total number of domains ',npni_opt*npnj_opt
+      WRITE(numout,*) ' '
+      WRITE(numout,*) ' jpni=',npni_opt ,' jpnj=',npnj_opt
+      WRITE(numout,*) ' jpi= ',jpimax_opt ,' jpj= ',jpjmax_opt
+      WRITE(numout,*) 
+      WRITE(numout,*) ' Number of ocean processors  ', npni_opt*npnj_opt-nland_opt
+      WRITE(numout,*) ' Number of land processors   ', nland_opt
+      WRITE(numout,*) ' Mean ocean coverage         ', float(ioce_opt)/float(npni_opt*npnj_opt-nland_opt)/float(jpimax_opt*jpjmax_opt)
+      WRITE(numout,*) ' Minimum ocean coverage      ', oce_cover_min_opt
+      WRITE(numout,*) ' Maximum ocean coverage      ', oce_cover_max_opt
+      WRITE(numout,*) ' nb of proc with coverage         < 10 % ', inf10_opt
+      WRITE(numout,*) ' nb of proc with coverage 10 < nb < 30 % ', inf30_opt - inf10_opt
+      WRITE(numout,*) ' nb of proc with coverage 30 < nb < 50 % ', inf50_opt - inf30_opt
+      WRITE(numout,*) ' Number of computed points   ', (npni_opt*npnj_opt-nland_opt)*jpimax_opt*jpjmax_opt
+      WRITE(numout,*) ' Overhead of computed points ', (npni_opt*npnj_opt-nland_opt)*jpimax_opt*jpjmax_opt-npiglo*npjglo
+      WRITE(numout,*) ' % sup (computed / global)   ', ratio_opt
+      WRITE(numout,*)
+   ENDIF
+   CLOSE(numout)
+   !
+   STOP
+END PROGRAM mpp_optimize
