Changeset 10330


Ignore:
Timestamp:
2018-11-19T10:41:29+01:00 (21 months ago)
Author:
smasson
Message:

dev_r10164_HPC09_ESIWACE_PREP_MERGE: action 6: find the best mpi domain decomposition, see #2133

Location:
NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/IOM/in_out_manager.F90

    r10068 r10330  
    122122   INTEGER ::   numtime         =   -1      !: logical unit for timing 
    123123   INTEGER ::   numout          =    6      !: logical unit for output print; Set to stdout to ensure any 
     124   INTEGER ::   numnul          =   -1      !: logical unit for /dev/null 
    124125      !                                     !  early output can be collected; do not change 
    125126   INTEGER ::   numnam_ref      =   -1      !: logical unit for reference namelist 
  • NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/LBC/lib_mpp.F90

    r10314 r10330  
    8686   PUBLIC   mppscatter, mppgather 
    8787   PUBLIC   mpp_ini_znl 
    88    PUBLIC   mppsize 
    8988   PUBLIC   mppsend, mpprecv                          ! needed by TAM and ICB routines 
    9089   PUBLIC   mpp_lnk_bdy_2d, mpp_lnk_bdy_3d, mpp_lnk_bdy_4d 
    91    PUBLIC   mpprank 
    9290    
    9391   !! * Interfaces 
     
    123121   INTEGER, PARAMETER         ::   nprocmax = 2**10   ! maximun dimension (required to be a power of 2) 
    124122 
    125    INTEGER ::   mppsize        ! number of process 
    126    INTEGER ::   mpprank        ! process number  [ 0 - size-1 ] 
     123   INTEGER, PUBLIC ::   mppsize        ! number of process 
     124   INTEGER, PUBLIC ::   mpprank        ! process number  [ 0 - size-1 ] 
    127125!$AGRIF_DO_NOT_TREAT 
    128126   INTEGER, PUBLIC ::   mpi_comm_oce   ! opa local communicator 
     
    220218      WRITE(ldtxt(ii),*) '      size exported buffer   nn_buffer   = ', nn_buffer,' bytes';   ii = ii + 1 
    221219      ! 
    222 #if defined key_agrif 
    223       IF( .NOT. Agrif_Root() ) THEN 
    224          jpni  = Agrif_Parent(jpni ) 
    225          jpnj  = Agrif_Parent(jpnj ) 
    226          jpnij = Agrif_Parent(jpnij) 
    227       ENDIF 
    228 #endif 
    229       ! 
    230       IF( jpnij < 1 ) THEN         ! If jpnij is not specified in namelist then we calculate it 
    231          jpnij = jpni * jpnj       ! this means there will be no land cutting out. 
    232       ENDIF 
    233  
    234220      IF( jpni < 1 .OR. jpnj < 1  ) THEN 
    235221         WRITE(ldtxt(ii),*) '      jpni, jpnj and jpnij will be calculated automatically' ;   ii = ii + 1 
     
    237223         WRITE(ldtxt(ii),*) '      processor grid extent in i         jpni = ',jpni       ;   ii = ii + 1 
    238224         WRITE(ldtxt(ii),*) '      processor grid extent in j         jpnj = ',jpnj       ;   ii = ii + 1 
    239          WRITE(ldtxt(ii),*) '      number of local domains           jpnij = ',jpnij      ;   ii = ii + 1 
    240225      ENDIF 
    241226 
     
    765750 
    766751 
    767    SUBROUTINE mppstop 
     752   SUBROUTINE mppstop( ldfinal )  
    768753      !!---------------------------------------------------------------------- 
    769754      !!                  ***  routine mppstop  *** 
     
    772757      !! 
    773758      !!---------------------------------------------------------------------- 
     759      LOGICAL, OPTIONAL, INTENT(in) :: ldfinal    ! source process number 
     760      LOGICAL ::   llfinal 
    774761      INTEGER ::   info 
    775762      !!---------------------------------------------------------------------- 
     
    777764      CALL mppsync 
    778765      CALL mpi_finalize( info ) 
     766      llfinal = .FALSE. 
     767      IF( PRESENT(ldfinal) ) llfinal = ldfinal 
     768      IF( .NOT. llfinal ) STOP 123456 
    779769      ! 
    780770   END SUBROUTINE mppstop 
     
    16821672      IF(lwp) THEN 
    16831673         WRITE(numout,cform_err) 
    1684          IF( PRESENT(cd1 ) )   WRITE(numout,*) cd1 
    1685          IF( PRESENT(cd2 ) )   WRITE(numout,*) cd2 
    1686          IF( PRESENT(cd3 ) )   WRITE(numout,*) cd3 
    1687          IF( PRESENT(cd4 ) )   WRITE(numout,*) cd4 
    1688          IF( PRESENT(cd5 ) )   WRITE(numout,*) cd5 
    1689          IF( PRESENT(cd6 ) )   WRITE(numout,*) cd6 
    1690          IF( PRESENT(cd7 ) )   WRITE(numout,*) cd7 
    1691          IF( PRESENT(cd8 ) )   WRITE(numout,*) cd8 
    1692          IF( PRESENT(cd9 ) )   WRITE(numout,*) cd9 
    1693          IF( PRESENT(cd10) )   WRITE(numout,*) cd10 
     1674         IF( PRESENT(cd1 ) )   WRITE(numout,*) TRIM(cd1) 
     1675         IF( PRESENT(cd2 ) )   WRITE(numout,*) TRIM(cd2) 
     1676         IF( PRESENT(cd3 ) )   WRITE(numout,*) TRIM(cd3) 
     1677         IF( PRESENT(cd4 ) )   WRITE(numout,*) TRIM(cd4) 
     1678         IF( PRESENT(cd5 ) )   WRITE(numout,*) TRIM(cd5) 
     1679         IF( PRESENT(cd6 ) )   WRITE(numout,*) TRIM(cd6) 
     1680         IF( PRESENT(cd7 ) )   WRITE(numout,*) TRIM(cd7) 
     1681         IF( PRESENT(cd8 ) )   WRITE(numout,*) TRIM(cd8) 
     1682         IF( PRESENT(cd9 ) )   WRITE(numout,*) TRIM(cd9) 
     1683         IF( PRESENT(cd10) )   WRITE(numout,*) TRIM(cd10) 
    16941684      ENDIF 
    16951685                               CALL FLUSH(numout    ) 
     
    17211711      IF(lwp) THEN 
    17221712         WRITE(numout,cform_war) 
    1723          IF( PRESENT(cd1 ) ) WRITE(numout,*) cd1 
    1724          IF( PRESENT(cd2 ) ) WRITE(numout,*) cd2 
    1725          IF( PRESENT(cd3 ) ) WRITE(numout,*) cd3 
    1726          IF( PRESENT(cd4 ) ) WRITE(numout,*) cd4 
    1727          IF( PRESENT(cd5 ) ) WRITE(numout,*) cd5 
    1728          IF( PRESENT(cd6 ) ) WRITE(numout,*) cd6 
    1729          IF( PRESENT(cd7 ) ) WRITE(numout,*) cd7 
    1730          IF( PRESENT(cd8 ) ) WRITE(numout,*) cd8 
    1731          IF( PRESENT(cd9 ) ) WRITE(numout,*) cd9 
    1732          IF( PRESENT(cd10) ) WRITE(numout,*) cd10 
     1713         IF( PRESENT(cd1 ) ) WRITE(numout,*) TRIM(cd1) 
     1714         IF( PRESENT(cd2 ) ) WRITE(numout,*) TRIM(cd2) 
     1715         IF( PRESENT(cd3 ) ) WRITE(numout,*) TRIM(cd3) 
     1716         IF( PRESENT(cd4 ) ) WRITE(numout,*) TRIM(cd4) 
     1717         IF( PRESENT(cd5 ) ) WRITE(numout,*) TRIM(cd5) 
     1718         IF( PRESENT(cd6 ) ) WRITE(numout,*) TRIM(cd6) 
     1719         IF( PRESENT(cd7 ) ) WRITE(numout,*) TRIM(cd7) 
     1720         IF( PRESENT(cd8 ) ) WRITE(numout,*) TRIM(cd8) 
     1721         IF( PRESENT(cd9 ) ) WRITE(numout,*) TRIM(cd9) 
     1722         IF( PRESENT(cd10) ) WRITE(numout,*) TRIM(cd10) 
    17331723      ENDIF 
    17341724      CALL FLUSH(numout) 
     
    17711761      knum=get_unit() 
    17721762#endif 
     1763      IF( TRIM(cdfile) == '/dev/null' )   clfile = TRIM(cdfile)   ! force the use of /dev/null 
    17731764      ! 
    17741765      iost=0 
     
    17781769         OPEN( UNIT=knum, FILE=clfile, FORM=cdform, ACCESS=cdacce, STATUS=cdstat             , ERR=100, IOSTAT=iost ) 
    17791770      ENDIF 
     1771      IF( iost /= 0 .AND. TRIM(clfile) == '/dev/null' ) & 
     1772         &  OPEN(UNIT=knum,FILE='NUL', FORM=cdform, ACCESS=cdacce, STATUS=cdstat             , ERR=100, IOSTAT=iost )   ! for windows 
    17801773      IF( iost == 0 ) THEN 
    17811774         IF(ldwp) THEN 
  • NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/LBC/mppini.F90

    r10297 r10330  
    3636   PUBLIC mpp_init       ! called by opa.F90 
    3737 
     38   INTEGER :: numbot = -1  ! 'bottom_level' local logical unit 
     39   INTEGER :: numbdy = -1  ! 'bdy_msk'      local logical unit 
     40    
    3841   !!---------------------------------------------------------------------- 
    3942   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    136139      !!---------------------------------------------------------------------- 
    137140      INTEGER ::   ji, jj, jn, jproc, jarea   ! dummy loop indices 
     141      INTEGER ::   inij, inijmin 
     142      INTEGER ::   i2add 
    138143      INTEGER ::   inum                       ! local logical unit 
    139       INTEGER ::   idir, ifreq, icont, isurf  ! local integers 
     144      INTEGER ::   idir, ifreq, icont         ! local integers 
    140145      INTEGER ::   ii, il1, ili, imil         !   -       - 
    141146      INTEGER ::   ij, il2, ilj, ijm1         !   -       - 
    142147      INTEGER ::   iino, ijno, iiso, ijso     !   -       - 
    143148      INTEGER ::   iiea, ijea, iiwe, ijwe     !   -       - 
    144       INTEGER ::   iresti, irestj, iarea0     !   -       - 
    145       INTEGER ::   ierr                       ! local logical unit 
    146       REAL(wp)::   zidom, zjdom               ! local scalars 
     149      INTEGER ::   iarea0                     !   -       - 
     150      INTEGER ::   ierr, ios                  !  
    147151      INTEGER, ALLOCATABLE, DIMENSION(:)     ::   iin, ii_nono, ii_noea          ! 1D workspace 
    148152      INTEGER, ALLOCATABLE, DIMENSION(:)     ::   ijn, ii_noso, ii_nowe          !  -     - 
     
    151155      INTEGER, ALLOCATABLE, DIMENSION(:,:) ::   ilei, ildi, iono, ioea         !  -     - 
    152156      INTEGER, ALLOCATABLE, DIMENSION(:,:) ::   ilej, ildj, ioso, iowe         !  -     - 
    153       INTEGER, DIMENSION(jpiglo,jpjglo) ::   imask   ! 2D global domain workspace 
    154       !!---------------------------------------------------------------------- 
    155  
     157      LOGICAL, ALLOCATABLE, DIMENSION(:,:) ::   llisoce                        !  -     - 
     158      NAMELIST/nambdy/ ln_bdy, nb_bdy, ln_coords_file, cn_coords_file,           & 
     159           &             ln_mask_file, cn_mask_file, cn_dyn2d, nn_dyn2d_dta,     & 
     160           &             cn_dyn3d, nn_dyn3d_dta, cn_tra, nn_tra_dta,             &   
     161           &             ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp, rn_time_dmp_out, & 
     162           &             cn_ice, nn_ice_dta,                                     & 
     163           &             rn_ice_tem, rn_ice_sal, rn_ice_age,                     & 
     164           &             ln_vol, nn_volctl, nn_rimwidth, nb_jpk_bdy 
     165      !!---------------------------------------------------------------------- 
     166 
     167      ! do we need to take into account bdy_msk? 
     168      REWIND( numnam_ref )              ! Namelist nambdy in reference namelist : BDY 
     169      READ  ( numnam_ref, nambdy, IOSTAT = ios, ERR = 903) 
     170903   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy in reference namelist (mppini)', lwp ) 
     171      REWIND( numnam_cfg )              ! Namelist nambdy in configuration namelist : BDY 
     172      READ  ( numnam_cfg, nambdy, IOSTAT = ios, ERR = 904 ) 
     173904   IF( ios >  0 )   CALL ctl_nam ( ios , 'nambdy in configuration namelist (mppini)', lwp ) 
     174      ! 
     175      IF(               ln_read_cfg ) CALL iom_open( cn_domcfg,    numbot ) 
     176      IF( ln_bdy .AND. ln_mask_file ) CALL iom_open( cn_mask_file, numbdy ) 
     177      ! 
     178      !  1. Dimension arrays for subdomains 
     179      ! ----------------------------------- 
     180      ! 
    156181      ! If dimensions of processor grid weren't specified in the namelist file 
    157182      ! then we calculate them here now that we have our communicator size 
    158       IF( jpni < 1 .OR. jpnj < 1 )   CALL mpp_init_partition( mppsize ) 
    159       ! 
    160 #if defined key_agrif 
    161       IF( jpnij /= jpni*jpnj ) CALL ctl_stop( 'STOP', 'Cannot remove land proc with AGRIF' ) 
    162 #endif 
    163       ! 
     183      IF( jpni < 1 .OR. jpnj < 1 )   CALL mpp_init_bestpartition( mppsize, jpni, jpnj ) 
     184       
     185      ! look for land mpi subdomains... 
     186      ALLOCATE( llisoce(jpni,jpnj) ) 
     187      CALL mpp_init_isoce( jpni, jpnj, llisoce ) 
     188      inijmin = COUNT( llisoce )   ! number of oce subdomains 
     189 
     190      IF( mppsize < inijmin ) THEN 
     191         WRITE(ctmp1,*) '   With this specified domain decomposition: jpni =', jpni, ' jpnj =', jpnj 
     192         WRITE(ctmp2,*) '   we can eliminate only ', jpni*jpnj - inijmin, ' land mpi subdomains therefore ' 
     193         WRITE(ctmp3,*) '   the number of ocean mpi subdomains (', inijmin,') exceed the number of MPI processes:', mppsize 
     194         WRITE(ctmp4,*) '   ==>>> There is the list of best domain decompositions you should use: ' 
     195         CALL ctl_stop( 'mpp_init:', '~~~~~~~~ ', ctmp1, ctmp2, ctmp3, ctmp4 ) 
     196         CALL mpp_init_bestpartition( mppsize, ldlist = .TRUE. )   ! must be done by all core 
     197         CALL ctl_stop( 'STOP' ) 
     198      ENDIF 
     199 
     200      IF( mppsize > jpni*jpnj ) THEN 
     201         WRITE(ctmp1,*) '   The number of mpi processes: ', mppsize 
     202         WRITE(ctmp2,*) '   exceeds the maximum number of subdomains (ocean+land) = ', jpni*jpnj 
     203         WRITE(ctmp3,*) '   defined by the following domain decomposition: jpni =', jpni, ' jpnj =', jpnj 
     204         WRITE(ctmp4,*) '   ==>>> There is the list of best domain decompositions you should use: ' 
     205         CALL ctl_stop( 'mpp_init:', '~~~~~~~~ ', ctmp1, ctmp2, ctmp3, ctmp4 ) 
     206         CALL mpp_init_bestpartition( mppsize, ldlist = .TRUE. )   ! must be done by all core 
     207         CALL ctl_stop( 'STOP' ) 
     208      ENDIF 
     209 
     210      jpnij = mppsize   ! force jpnij definition <-- remove as much land subdomains as needed to reach this condition 
     211      IF( mppsize > inijmin ) THEN 
     212         WRITE(ctmp1,*) '   The number of mpi processes: ', mppsize 
     213         WRITE(ctmp2,*) '   exceeds the maximum number of ocean subdomains = ', inijmin 
     214         WRITE(ctmp3,*) '   we suppressed ', jpni*jpnj - mppsize, ' land subdomains ' 
     215         WRITE(ctmp4,*) '   BUT we had to keep ', mppsize - inijmin, ' land subdomains that are useless...' 
     216         CALL ctl_warn( 'mpp_init:', '~~~~~~~~ ', ctmp1, ctmp2, ctmp3, ctmp4 ) 
     217      ELSE   ! mppsize = inijmin 
     218         IF(lwp) THEN 
     219            WRITE(numout,*) 'mpp_init: You use an optimal domain decomposition' 
     220            WRITE(numout,*) '~~~~~~~~ ' 
     221            WRITE(numout,*) '   Number of mpi processes: ', mppsize 
     222            WRITE(numout,*) '   Number of ocean subdomains = ', inijmin 
     223            WRITE(numout,*) '   Number of suppressed land subdomains = ', jpni*jpnj - inijmin 
     224         ENDIF 
     225      ENDIF 
     226 
     227      IF( numbot /= -1 )   CALL iom_close( numbot ) 
     228      IF( numbdy /= -1 )   CALL iom_close( numbdy ) 
     229     
    164230      ALLOCATE(  nfiimpp(jpni,jpnj), nfipproc(jpni,jpnj), nfilcit(jpni,jpnj) ,    & 
    165231         &       nimppt(jpnij) , ibonit(jpnij) , nlcit(jpnij) , nlcjt(jpnij) ,    & 
     
    176242      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'mpp_init: unable to allocate standard ocean arrays' ) 
    177243       
    178       ! 
    179244#if defined key_agrif 
    180245      IF( .NOT. Agrif_Root() ) THEN       ! AGRIF children: specific setting (cf. agrif_user.F90) 
     
    186251      ENDIF 
    187252#endif 
    188  
    189 #if defined key_nemocice_decomp 
    190       jpimax = ( nx_global+2-2*nn_hls + (jpni-1) ) / jpni + 2*nn_hls    ! first  dim. 
    191       jpjmax = ( ny_global+2-2*nn_hls + (jpnj-1) ) / jpnj + 2*nn_hls    ! second dim.  
    192 #else 
    193       jpimax = ( jpiglo - 2*nn_hls + (jpni-1) ) / jpni + 2*nn_hls    ! first  dim. 
    194       jpjmax = ( jpjglo - 2*nn_hls + (jpnj-1) ) / jpnj + 2*nn_hls    ! second dim. 
    195 #endif 
    196  
    197       ! 
    198       IF ( jpni * jpnj == jpnij ) THEN    ! regular domain lay out over processors 
    199          imask(:,:) = 1                
    200       ELSEIF ( jpni*jpnj > jpnij ) THEN   ! remove land-only processor (i.e. where imask(:,:)=0) 
    201          CALL mpp_init_mask( imask )    
    202       ELSE                                ! error 
    203          CALL ctl_stop( 'mpp_init: jpnij > jpni x jpnj. Check namelist setting!' ) 
    204       ENDIF 
    205       ! 
    206       !  1. Dimension arrays for subdomains 
     253      ! 
     254      !  2. Index arrays for subdomains 
    207255      ! ----------------------------------- 
    208       !  Computation of local domain sizes ilci() ilcj() 
    209       !  These dimensions depend on global sizes jpni,jpnj and jpiglo,jpjglo 
    210       !  The subdomains are squares lesser than or equal to the global 
    211       !  dimensions divided by the number of processors minus the overlap array. 
    212256      ! 
    213257      nreci = 2 * nn_hls 
    214258      nrecj = 2 * nn_hls 
    215       iresti = 1 + MOD( jpiglo - nreci -1 , jpni ) 
    216       irestj = 1 + MOD( jpjglo - nrecj -1 , jpnj ) 
    217       ! 
    218       !  Need to use jpimax and jpjmax here since jpi and jpj not yet defined 
    219 #if defined key_nemocice_decomp 
    220       ! Change padding to be consistent with CICE 
    221       ilci(1:jpni-1      ,:) = jpimax 
    222       ilci(jpni          ,:) = jpiglo - (jpni - 1) * (jpimax - nreci) 
    223       ! 
    224       ilcj(:,      1:jpnj-1) = jpjmax 
    225       ilcj(:,          jpnj) = jpjglo - (jpnj - 1) * (jpjmax - nrecj) 
    226 #else 
    227       ilci(1:iresti      ,:) = jpimax 
    228       ilci(iresti+1:jpni ,:) = jpimax-1 
    229  
    230       ilcj(:,      1:irestj) = jpjmax 
    231       ilcj(:, irestj+1:jpnj) = jpjmax-1 
    232 #endif 
    233       ! 
    234       zidom = nreci + sum(ilci(:,1) - nreci ) 
    235       zjdom = nrecj + sum(ilcj(1,:) - nrecj ) 
     259      CALL mpp_basic_decomposition( jpni, jpnj, jpimax, jpjmax, iimppt, ijmppt, ilci, ilcj ) 
     260      nfiimpp(:,:) = iimppt(:,:) 
     261      nfilcit(:,:) = ilci(:,:) 
    236262      ! 
    237263      IF(lwp) THEN 
    238264         WRITE(numout,*) 
    239          WRITE(numout,*) 'mpp_init : MPI Message Passing MPI - domain lay out over processors' 
    240          WRITE(numout,*) '~~~~~~~~ ' 
     265         WRITE(numout,*) 'MPI Message Passing MPI - domain lay out over processors' 
     266         WRITE(numout,*) 
    241267         WRITE(numout,*) '   defines mpp subdomains' 
    242          WRITE(numout,*) '      iresti = ', iresti, ' jpni = ', jpni   
    243          WRITE(numout,*) '      irestj = ', irestj, ' jpnj = ', jpnj 
     268         WRITE(numout,*) '      jpni = ', jpni   
     269         WRITE(numout,*) '      jpnj = ', jpnj 
    244270         WRITE(numout,*) 
    245          WRITE(numout,*) '      sum ilci(i,1) = ', zidom, ' jpiglo = ', jpiglo 
    246          WRITE(numout,*) '      sum ilcj(1,j) = ', zjdom, ' jpjglo = ', jpjglo 
    247       ENDIF 
    248  
    249       !  2. Index arrays for subdomains 
    250       ! ------------------------------- 
    251       iimppt(:,:) =  1 
    252       ijmppt(:,:) =  1 
    253       ipproc(:,:) = -1 
    254       ! 
    255       IF( jpni > 1 ) THEN 
    256          DO jj = 1, jpnj 
    257             DO ji = 2, jpni 
    258                iimppt(ji,jj) = iimppt(ji-1,jj) + ilci(ji-1,jj) - nreci 
    259             END DO 
    260          END DO 
    261       ENDIF 
    262       nfiimpp(:,:) = iimppt(:,:) 
    263       ! 
    264       IF( jpnj > 1 )THEN 
    265          DO jj = 2, jpnj 
    266             DO ji = 1, jpni 
    267                ijmppt(ji,jj) = ijmppt(ji,jj-1) + ilcj(ji,jj-1) - nrecj 
    268             END DO 
    269          END DO 
    270       ENDIF 
    271  
     271         WRITE(numout,*) '      sum ilci(i,1) = ', sum(ilci(:,1)), ' jpiglo = ', jpiglo 
     272         WRITE(numout,*) '      sum ilcj(1,j) = ', sum(ilcj(1,:)), ' jpjglo = ', jpjglo 
     273      ENDIF 
     274      
    272275      ! 3. Subdomain description in the Regular Case 
    273276      ! -------------------------------------------- 
     
    277280      l_Jperio = jpnj == 1 .AND. (jperio == 2 .OR. jperio == 7) 
    278281       
    279       icont = -1 
    280282      DO jarea = 1, jpni*jpnj 
     283         ! 
    281284         iarea0 = jarea - 1 
    282285         ii = 1 + MOD(iarea0,jpni) 
     
    334337         ENDIF 
    335338         ! 
    336          ! Check wet points over the entire domain to preserve the MPI communication stencil 
    337          isurf = 0 
    338          DO jj = 1, ilj 
    339             DO  ji = 1, ili 
    340                IF( imask(ji+iimppt(ii,ij)-1, jj+ijmppt(ii,ij)-1) == 1)   isurf = isurf+1 
    341             END DO 
    342          END DO 
    343          ! 
    344          IF( isurf /= 0 ) THEN 
     339      END DO 
     340 
     341      ! 4. deal with land subdomains 
     342      ! ---------------------------- 
     343      ! 
     344      ! specify which subdomains are oce subdomains; other are land subdomains 
     345      ipproc(:,:) = -1 
     346      icont = -1 
     347      DO jarea = 1, jpni*jpnj 
     348         iarea0 = jarea - 1 
     349         ii = 1 + MOD(iarea0,jpni) 
     350         ij = 1 +     iarea0/jpni 
     351         IF( llisoce(ii,ij) ) THEN 
    345352            icont = icont + 1 
    346353            ipproc(ii,ij) = icont 
     
    349356         ENDIF 
    350357      END DO 
    351       ! 
     358      ! if needed add some land subdomains to reach jpnij active subdomains 
     359      i2add = jpnij - inijmin 
     360      DO jarea = 1, jpni*jpnj 
     361         iarea0 = jarea - 1 
     362         ii = 1 + MOD(iarea0,jpni) 
     363         ij = 1 +     iarea0/jpni 
     364         IF( .NOT. llisoce(ii,ij) .AND. i2add > 0 ) THEN 
     365            icont = icont + 1 
     366            ipproc(ii,ij) = icont 
     367            iin(icont+1) = ii 
     368            ijn(icont+1) = ij 
     369            i2add = i2add - 1 
     370         ENDIF 
     371      END DO 
    352372      nfipproc(:,:) = ipproc(:,:) 
    353373 
    354       ! Check potential error 
    355       IF( icont+1 /= jpnij ) THEN 
    356          WRITE(ctmp1,*) ' jpni =',jpni,' jpnj =',jpnj 
    357          WRITE(ctmp2,*) ' jpnij =',jpnij, '< jpni x jpnj'  
    358          WRITE(ctmp3,*) ' ***********, mpp_init2 finds jpnij=',icont+1 
    359          CALL ctl_stop( 'STOP', 'mpp_init: Eliminate land processors algorithm', '', ctmp1, ctmp2, '', ctmp3 ) 
    360       ENDIF 
    361  
    362       ! 4. Subdomain print 
     374      ! neighbour treatment: change ibondi, ibondj if next to a land zone 
     375      DO jarea = 1, jpni*jpnj 
     376         ii = 1 + MOD( jarea-1  , jpni ) 
     377         ij = 1 +     (jarea-1) / jpni 
     378         ! land-only area with an active n neigbour 
     379         IF ( ipproc(ii,ij) == -1 .AND. 0 <= iono(ii,ij) .AND. iono(ii,ij) <= jpni*jpnj-1 ) THEN 
     380            iino = 1 + MOD( iono(ii,ij) , jpni )                    ! ii index of this n neigbour 
     381            ijno = 1 +      iono(ii,ij) / jpni                      ! ij index of this n neigbour 
     382            ! In case of north fold exchange: I am the n neigbour of my n neigbour!! (#1057) 
     383            ! --> for northern neighbours of northern row processors (in case of north-fold) 
     384            !     need to reverse the LOGICAL direction of communication  
     385            idir = 1                                           ! we are indeed the s neigbour of this n neigbour 
     386            IF( ij == jpnj .AND. ijno == jpnj )   idir = -1    ! both are on the last row, we are in fact the n neigbour 
     387            IF( ibondj(iino,ijno) == idir     )   ibondj(iino,ijno) =   2     ! this n neigbour had only a s/n neigbour -> no more 
     388            IF( ibondj(iino,ijno) == 0        )   ibondj(iino,ijno) = -idir   ! this n neigbour had both, n-s neighbours -> keep 1 
     389         ENDIF 
     390         ! land-only area with an active s neigbour 
     391         IF( ipproc(ii,ij) == -1 .AND. 0 <= ioso(ii,ij) .AND. ioso(ii,ij) <= jpni*jpnj-1 ) THEN 
     392            iiso = 1 + MOD( ioso(ii,ij) , jpni )                    ! ii index of this s neigbour 
     393            ijso = 1 +      ioso(ii,ij) / jpni                      ! ij index of this s neigbour 
     394            IF( ibondj(iiso,ijso) == -1 )   ibondj(iiso,ijso) = 2   ! this s neigbour had only a n neigbour    -> no more neigbour 
     395            IF( ibondj(iiso,ijso) ==  0 )   ibondj(iiso,ijso) = 1   ! this s neigbour had both, n-s neighbours -> keep s neigbour 
     396         ENDIF 
     397         ! land-only area with an active e neigbour 
     398         IF( ipproc(ii,ij) == -1 .AND. 0 <= ioea(ii,ij) .AND. ioea(ii,ij) <= jpni*jpnj-1 ) THEN 
     399            iiea = 1 + MOD( ioea(ii,ij) , jpni )                    ! ii index of this e neigbour 
     400            ijea = 1 +      ioea(ii,ij) / jpni                      ! ij index of this e neigbour 
     401            IF( ibondi(iiea,ijea) == 1 )   ibondi(iiea,ijea) =  2   ! this e neigbour had only a w neigbour    -> no more neigbour 
     402            IF( ibondi(iiea,ijea) == 0 )   ibondi(iiea,ijea) = -1   ! this e neigbour had both, e-w neighbours -> keep e neigbour 
     403         ENDIF 
     404         ! land-only area with an active w neigbour 
     405         IF( ipproc(ii,ij) == -1 .AND. 0 <= iowe(ii,ij) .AND. iowe(ii,ij) <= jpni*jpnj-1) THEN 
     406            iiwe = 1 + MOD( iowe(ii,ij) , jpni )                    ! ii index of this w neigbour 
     407            ijwe = 1 +      iowe(ii,ij) / jpni                      ! ij index of this w neigbour 
     408            IF( ibondi(iiwe,ijwe) == -1 )   ibondi(iiwe,ijwe) = 2   ! this w neigbour had only a e neigbour    -> no more neigbour 
     409            IF( ibondi(iiwe,ijwe) ==  0 )   ibondi(iiwe,ijwe) = 1   ! this w neigbour had both, e-w neighbours -> keep w neigbour 
     410         ENDIF 
     411      END DO 
     412 
     413      ! Update il[de][ij] according to modified ibond[ij] 
     414      ! ---------------------- 
     415      DO jproc = 1, jpnij 
     416         ii = iin(jproc) 
     417         ij = ijn(jproc) 
     418         IF( ibondi(ii,ij) == -1 .OR. ibondi(ii,ij) == 2 ) ildi(ii,ij) =  1 
     419         IF( ibondi(ii,ij) ==  1 .OR. ibondi(ii,ij) == 2 ) ilei(ii,ij) = ilci(ii,ij) 
     420         IF( ibondj(ii,ij) == -1 .OR. ibondj(ii,ij) == 2 ) ildj(ii,ij) =  1 
     421         IF( ibondj(ii,ij) ==  1 .OR. ibondj(ii,ij) == 2 ) ilej(ii,ij) = ilcj(ii,ij) 
     422      END DO 
     423       
     424      ! 5. Subdomain print 
    363425      ! ------------------ 
    364426      IF(lwp) THEN 
     
    385447 9404    FORMAT('           *  '   ,20('      ',i3,'   *   ') ) 
    386448      ENDIF 
    387  
    388       ! 5. neighbour treatment: change ibondi, ibondj if next to a land zone 
    389       ! ---------------------- 
    390       DO jarea = 1, jpni*jpnj 
    391          ii = 1 + MOD( jarea-1  , jpni ) 
    392          ij = 1 +     (jarea-1) / jpni 
    393          ! land-only area with an active n neigbour 
    394          IF ( ipproc(ii,ij) == -1 .AND. 0 <= iono(ii,ij) .AND. iono(ii,ij) <= jpni*jpnj-1 ) THEN 
    395             iino = 1 + MOD( iono(ii,ij) , jpni )                    ! ii index of this n neigbour 
    396             ijno = 1 +      iono(ii,ij) / jpni                      ! ij index of this n neigbour 
    397             ! In case of north fold exchange: I am the n neigbour of my n neigbour!! (#1057) 
    398             ! --> for northern neighbours of northern row processors (in case of north-fold) 
    399             !     need to reverse the LOGICAL direction of communication  
    400             idir = 1                                           ! we are indeed the s neigbour of this n neigbour 
    401             IF( ij == jpnj .AND. ijno == jpnj )   idir = -1    ! both are on the last row, we are in fact the n neigbour 
    402             IF( ibondj(iino,ijno) == idir     )   ibondj(iino,ijno) =   2     ! this n neigbour had only a s/n neigbour -> no more 
    403             IF( ibondj(iino,ijno) == 0        )   ibondj(iino,ijno) = -idir   ! this n neigbour had both, n-s neighbours -> keep 1 
    404          ENDIF 
    405          ! land-only area with an active s neigbour 
    406          IF( ipproc(ii,ij) == -1 .AND. 0 <= ioso(ii,ij) .AND. ioso(ii,ij) <= jpni*jpnj-1 ) THEN 
    407             iiso = 1 + MOD( ioso(ii,ij) , jpni )                    ! ii index of this s neigbour 
    408             ijso = 1 +      ioso(ii,ij) / jpni                      ! ij index of this s neigbour 
    409             IF( ibondj(iiso,ijso) == -1 )   ibondj(iiso,ijso) = 2   ! this s neigbour had only a n neigbour    -> no more neigbour 
    410             IF( ibondj(iiso,ijso) ==  0 )   ibondj(iiso,ijso) = 1   ! this s neigbour had both, n-s neighbours -> keep s neigbour 
    411          ENDIF 
    412          ! land-only area with an active e neigbour 
    413          IF( ipproc(ii,ij) == -1 .AND. 0 <= ioea(ii,ij) .AND. ioea(ii,ij) <= jpni*jpnj-1 ) THEN 
    414             iiea = 1 + MOD( ioea(ii,ij) , jpni )                    ! ii index of this e neigbour 
    415             ijea = 1 +      ioea(ii,ij) / jpni                      ! ij index of this e neigbour 
    416             IF( ibondi(iiea,ijea) == 1 )   ibondi(iiea,ijea) =  2   ! this e neigbour had only a w neigbour    -> no more neigbour 
    417             IF( ibondi(iiea,ijea) == 0 )   ibondi(iiea,ijea) = -1   ! this e neigbour had both, e-w neighbours -> keep e neigbour 
    418          ENDIF 
    419          ! land-only area with an active w neigbour 
    420          IF( ipproc(ii,ij) == -1 .AND. 0 <= iowe(ii,ij) .AND. iowe(ii,ij) <= jpni*jpnj-1) THEN 
    421             iiwe = 1 + MOD( iowe(ii,ij) , jpni )                    ! ii index of this w neigbour 
    422             ijwe = 1 +      iowe(ii,ij) / jpni                      ! ij index of this w neigbour 
    423             IF( ibondi(iiwe,ijwe) == -1 )   ibondi(iiwe,ijwe) = 2   ! this w neigbour had only a e neigbour    -> no more neigbour 
    424             IF( ibondi(iiwe,ijwe) ==  0 )   ibondi(iiwe,ijwe) = 1   ! this w neigbour had both, e-w neighbours -> keep w neigbour 
    425          ENDIF 
    426       END DO 
    427  
    428       ! Update il[de][ij] according to modified ibond[ij] 
    429       ! ---------------------- 
    430       DO jproc = 1, jpnij 
    431          ii = iin(jproc) 
    432          ij = ijn(jproc) 
    433          IF( ibondi(ii,ij) == -1 .OR. ibondi(ii,ij) == 2 ) ildi(ii,ij) =  1 
    434          IF( ibondi(ii,ij) ==  1 .OR. ibondi(ii,ij) == 2 ) ilei(ii,ij) = ilci(ii,ij) 
    435          IF( ibondj(ii,ij) == -1 .OR. ibondj(ii,ij) == 2 ) ildj(ii,ij) =  1 
    436          IF( ibondj(ii,ij) ==  1 .OR. ibondj(ii,ij) == 2 ) ilej(ii,ij) = ilcj(ii,ij) 
    437       END DO 
    438449          
    439450      ! just to save nono etc for all proc 
     
    516527         njmppt(jproc) = ijmppt(ii,ij)  
    517528      END DO 
    518       nfilcit(:,:) = ilci(:,:) 
    519529 
    520530      ! Save processor layout in ascii file 
     
    610620         &       iimppt, ijmppt, ibondi, ibondj, ipproc, ipolj,   & 
    611621         &       ilci, ilcj, ilei, ilej, ildi, ildj,              & 
    612          &       iono, ioea, ioso, iowe) 
     622         &       iono, ioea, ioso, iowe, llisoce) 
    613623      ! 
    614624    END SUBROUTINE mpp_init 
    615625 
    616626 
    617     SUBROUTINE mpp_init_mask( kmask ) 
    618       !!---------------------------------------------------------------------- 
    619       !!                  ***  ROUTINE mpp_init_mask  *** 
    620       !! 
    621       !! ** Purpose : Read relevant bathymetric information in a global array 
    622       !!              in order to provide a land/sea mask used for the elimination 
     627    SUBROUTINE mpp_basic_decomposition( knbi, knbj, kimax, kjmax, kimppt, kjmppt, klci, klcj) 
     628      !!---------------------------------------------------------------------- 
     629      !!                  ***  ROUTINE mpp_basic_decomposition  *** 
     630      !!                     
     631      !! ** Purpose :   Lay out the global domain over processors. 
     632      !! 
     633      !! ** Method  :   Global domain is distributed in smaller local domains. 
     634      !! 
     635      !! ** Action : - set for all knbi*knbj domains: 
     636      !!                    kimppt     : longitudinal index 
     637      !!                    kjmppt     : latitudinal  index 
     638      !!                    klci       : first dimension 
     639      !!                    klcj       : second dimension 
     640      !!---------------------------------------------------------------------- 
     641      INTEGER,                       INTENT(in   ) ::   knbi, knbj 
     642      INTEGER,                       INTENT(  out) ::   kimax, kjmax 
     643      INTEGER, DIMENSION(knbi,knbj), INTENT(  out) ::   kimppt, kjmppt 
     644      INTEGER, DIMENSION(knbi,knbj), INTENT(  out) ::   klci, klcj 
     645      ! 
     646      INTEGER ::   ji, jj 
     647      INTEGER ::   iresti, irestj 
     648      INTEGER ::   ireci, irecj 
     649      !!---------------------------------------------------------------------- 
     650      ! 
     651#if defined key_nemocice_decomp 
     652      jpimax = ( nx_global+2-2*nn_hls + (knbi-1) ) / knbi + 2*nn_hls    ! first  dim. 
     653      jpjmax = ( ny_global+2-2*nn_hls + (knbj-1) ) / knbj + 2*nn_hls    ! second dim.  
     654#else 
     655      kimax = ( jpiglo - 2*nn_hls + (knbi-1) ) / knbi + 2*nn_hls    ! first  dim. 
     656      kjmax = ( jpjglo - 2*nn_hls + (knbj-1) ) / knbj + 2*nn_hls    ! second dim. 
     657#endif 
     658      ! 
     659      !  1. Dimension arrays for subdomains 
     660      ! ----------------------------------- 
     661      !  Computation of local domain sizes klci() klcj() 
     662      !  These dimensions depend on global sizes knbi,knbj and jpiglo,jpjglo 
     663      !  The subdomains are squares lesser than or equal to the global 
     664      !  dimensions divided by the number of processors minus the overlap array. 
     665      ! 
     666      ireci = 2 * nn_hls 
     667      irecj = 2 * nn_hls 
     668      iresti = 1 + MOD( jpiglo - ireci -1 , knbi ) 
     669      irestj = 1 + MOD( jpjglo - irecj -1 , knbj ) 
     670      ! 
     671      !  Need to use kimax and kjmax here since jpi and jpj not yet defined 
     672#if defined key_nemocice_decomp 
     673      ! Change padding to be consistent with CICE 
     674      klci(1:knbi-1      ,:) = kimax 
     675      klci(knbi          ,:) = jpiglo - (knbi - 1) * (kimax - nreci) 
     676      klcj(:,      1:knbj-1) = kjmax 
     677      klcj(:,          knbj) = jpjglo - (knbj - 1) * (kjmax - nrecj) 
     678#else 
     679      klci(1:iresti      ,:) = kimax 
     680      klci(iresti+1:knbi ,:) = kimax-1 
     681      klcj(:,      1:irestj) = kjmax 
     682      klcj(:, irestj+1:knbj) = kjmax-1 
     683#endif 
     684 
     685      !  2. Index arrays for subdomains 
     686      ! ------------------------------- 
     687      kimppt(:,:) = 1 
     688      kjmppt(:,:) = 1 
     689      ! 
     690      IF( knbi > 1 ) THEN 
     691         DO jj = 1, knbj 
     692            DO ji = 2, knbi 
     693               kimppt(ji,jj) = kimppt(ji-1,jj) + klci(ji-1,jj) - ireci 
     694            END DO 
     695         END DO 
     696      ENDIF 
     697      ! 
     698      IF( knbj > 1 )THEN 
     699         DO jj = 2, knbj 
     700            DO ji = 1, knbi 
     701               kjmppt(ji,jj) = kjmppt(ji,jj-1) + klcj(ji,jj-1) - irecj 
     702            END DO 
     703         END DO 
     704      ENDIF 
     705       
     706   END SUBROUTINE mpp_basic_decomposition 
     707 
     708 
     709   SUBROUTINE mpp_init_bestpartition( knbij, knbi, knbj, ldlist ) 
     710      !!---------------------------------------------------------------------- 
     711      !!                 ***  ROUTINE mpp_init_bestpartition  *** 
     712      !! 
     713      !! ** Purpose : 
     714      !! 
     715      !! ** Method  : 
     716      !!---------------------------------------------------------------------- 
     717      INTEGER,           INTENT(in   ) ::   knbij         ! total number if subdomains               (knbi*knbj) 
     718      INTEGER, OPTIONAL, INTENT(  out) ::   knbi, knbj    ! number if subdomains along i and j (knbi and knbj) 
     719      LOGICAL, OPTIONAL, INTENT(in   ) ::   ldlist        ! .true.: print the list the best domain decompositions (with land) 
     720      ! 
     721      INTEGER :: ji, jj, ii, iitarget 
     722      INTEGER :: iszitst, iszjtst 
     723      INTEGER :: isziref, iszjref 
     724      INTEGER :: inbij, iszij 
     725      INTEGER :: inbimax, inbjmax, inbijmax 
     726      INTEGER :: isz0, isz1 
     727      INTEGER, DIMENSION(  :), ALLOCATABLE :: indexok 
     728      INTEGER, DIMENSION(  :), ALLOCATABLE :: inbi0, inbj0, inbij0   ! number of subdomains along i,j 
     729      INTEGER, DIMENSION(  :), ALLOCATABLE :: iszi0, iszj0, iszij0   ! max size of the subdomains along i,j 
     730      INTEGER, DIMENSION(  :), ALLOCATABLE :: inbi1, inbj1, inbij1   ! number of subdomains along i,j 
     731      INTEGER, DIMENSION(  :), ALLOCATABLE :: iszi1, iszj1, iszij1   ! max size of the subdomains along i,j 
     732      LOGICAL :: llist 
     733      LOGICAL, DIMENSION(:,:), ALLOCATABLE :: llmsk2d                 ! max size of the subdomains along i,j 
     734      LOGICAL, DIMENSION(:,:), ALLOCATABLE :: llisoce              !  -     - 
     735      REAL(wp)::   zpropland 
     736      !!---------------------------------------------------------------------- 
     737      ! 
     738      llist = .FALSE. 
     739      IF( PRESENT(ldlist) ) llist = ldlist 
     740 
     741      CALL mpp_init_landprop( zpropland )                      ! get the proportion of land point over the gloal domain 
     742      inbij = NINT( REAL(knbij, wp) / ( 1.0 - zpropland ) )    ! define the largest possible value for jpni*jpnj 
     743      ! 
     744      IF( llist ) THEN   ;   inbijmax = inbij*2 
     745      ELSE               ;   inbijmax = inbij 
     746      ENDIF 
     747      ! 
     748      ALLOCATE(inbi0(inbijmax),inbj0(inbijmax),iszi0(inbijmax),iszj0(inbijmax)) 
     749      ! 
     750      inbimax = 0 
     751      inbjmax = 0 
     752      isziref = jpiglo*jpjglo+1 
     753      iszjref = jpiglo*jpjglo+1 
     754      ! 
     755      ! get the list of knbi that gives a smaller jpimax than knbi-1 
     756      ! get the list of knbj that gives a smaller jpjmax than knbj-1 
     757      DO ji = 1, inbijmax       
     758#if defined key_nemocice_decomp 
     759         iszitst = ( nx_global+2-2*nn_hls + (ji-1) ) / ji + 2*nn_hls    ! first  dim. 
     760#else 
     761         iszitst = ( jpiglo - 2*nn_hls + (ji-1) ) / ji + 2*nn_hls 
     762#endif 
     763         IF( iszitst < isziref ) THEN 
     764            isziref = iszitst 
     765            inbimax = inbimax + 1 
     766            inbi0(inbimax) = ji 
     767            iszi0(inbimax) = isziref 
     768         ENDIF 
     769#if defined key_nemocice_decomp 
     770         iszjtst = ( ny_global+2-2*nn_hls + (ji-1) ) / ji + 2*nn_hls    ! first  dim. 
     771#else 
     772         iszjtst = ( jpjglo - 2*nn_hls + (ji-1) ) / ji + 2*nn_hls 
     773#endif 
     774         IF( iszjtst < iszjref ) THEN 
     775            iszjref = iszjtst 
     776            inbjmax = inbjmax + 1 
     777            inbj0(inbjmax) = ji 
     778            iszj0(inbjmax) = iszjref 
     779         ENDIF 
     780      END DO 
     781 
     782      ! combine these 2 lists to get all possible knbi*knbj <  inbijmax 
     783      ALLOCATE( llmsk2d(inbimax,inbjmax) ) 
     784      DO jj = 1, inbjmax 
     785         DO ji = 1, inbimax 
     786            IF ( inbi0(ji) * inbj0(jj) <= inbijmax ) THEN   ;   llmsk2d(ji,jj) = .TRUE. 
     787            ELSE                                            ;   llmsk2d(ji,jj) = .FALSE. 
     788            ENDIF 
     789         END DO 
     790      END DO 
     791      isz1 = COUNT(llmsk2d) 
     792      ALLOCATE( inbi1(isz1), inbj1(isz1), iszi1(isz1), iszj1(isz1) ) 
     793      ii = 0 
     794      DO jj = 1, inbjmax 
     795         DO ji = 1, inbimax 
     796            IF( llmsk2d(ji,jj) == .TRUE. ) THEN 
     797               ii = ii + 1 
     798               inbi1(ii) = inbi0(ji) 
     799               inbj1(ii) = inbj0(jj) 
     800               iszi1(ii) = iszi0(ji) 
     801               iszj1(ii) = iszj0(jj) 
     802            END IF 
     803         END DO 
     804      END DO 
     805      DEALLOCATE( inbi0, inbj0, iszi0, iszj0 ) 
     806      DEALLOCATE( llmsk2d ) 
     807 
     808      ALLOCATE( inbij1(isz1), iszij1(isz1) ) 
     809      inbij1(:) = inbi1(:) * inbj1(:) 
     810      iszij1(:) = iszi1(:) * iszj1(:) 
     811 
     812      ! if therr is no land and no print 
     813      IF( .NOT. llist .AND. numbot == -1 .AND. numbdy == -1 ) THEN 
     814         ! get the smaller partition which gives the smallest subdomain size 
     815         ii = MINLOC(inbij1, mask = iszij1 == MINVAL(iszij1), dim = 1) 
     816         knbi = inbi1(ii) 
     817         knbj = inbj1(ii) 
     818         DEALLOCATE( inbi1, inbj1, inbij1, iszi1, iszj1, iszij1 ) 
     819         RETURN 
     820      ENDIF 
     821 
     822      ! extract only the partitions which reduce the subdomain size in comparison with smaller partitions 
     823      ALLOCATE( indexok(isz1) )                                 ! to store indices of the best partitions 
     824      isz0 = 0                                                  ! number of best partitions      
     825      inbij = 1                                                 ! start with the min value of inbij1 => 1 
     826      iszij = jpiglo*jpjglo+1                                   ! default: larger than global domain 
     827      DO WHILE( inbij <= inbijmax )                             ! if we did not reach the max of inbij1 
     828         ii = MINLOC(iszij1, mask = inbij1 == inbij, dim = 1)   ! warning: send back the first occurence if multiple results 
     829         IF ( iszij1(ii) < iszij ) THEN 
     830            isz0 = isz0 + 1 
     831            indexok(isz0) = ii 
     832            iszij = iszij1(ii) 
     833         ENDIF 
     834         inbij = MINVAL(inbij1, mask = inbij1 > inbij)   ! warning: return largest integer value if mask = .false. everywhere 
     835      END DO 
     836      DEALLOCATE( inbij1, iszij1 ) 
     837 
     838      ! keep only the best partitions (sorted by increasing order of subdomains number and decreassing subdomain size) 
     839      ALLOCATE( inbi0(isz0), inbj0(isz0), iszi0(isz0), iszj0(isz0) ) 
     840      DO ji = 1, isz0 
     841         ii = indexok(ji) 
     842         inbi0(ji) = inbi1(ii) 
     843         inbj0(ji) = inbj1(ii) 
     844         iszi0(ji) = iszi1(ii) 
     845         iszj0(ji) = iszj1(ii) 
     846      END DO 
     847      DEALLOCATE( indexok, inbi1, inbj1, iszi1, iszj1 ) 
     848 
     849      IF( llist ) THEN  ! we print about 21 best partitions 
     850         IF(lwp) THEN 
     851            WRITE(numout,*) 
     852            WRITE(numout,         *) '                  For your information:' 
     853            WRITE(numout,'(a,i5,a)') '  list of the best partitions around ',   knbij, ' mpi processes' 
     854            WRITE(numout,         *) '  --------------------------------------', '-----', '--------------' 
     855            WRITE(numout,*) 
     856         END IF 
     857         iitarget = MINLOC( inbi0(:)*inbj0(:), mask = inbi0(:)*inbj0(:) >= knbij, dim = 1 ) 
     858         DO ji = MAX(1,iitarget-10), MIN(isz0,iitarget+10) 
     859            ALLOCATE( llisoce(inbi0(ji), inbj0(ji)) ) 
     860            CALL mpp_init_isoce( inbi0(ji), inbj0(ji), llisoce ) ! Warning: must be call by all cores (call mpp_sum) 
     861            inbij = COUNT(llisoce) 
     862            DEALLOCATE( llisoce ) 
     863            IF(lwp) WRITE(numout,'(a, i5, a, i5, a, i4, a, i4, a, i9, a, i5, a, i5, a)')    & 
     864               &     'nb_cores ' , inbij,' oce + ', inbi0(ji)*inbj0(ji) - inbij             & 
     865               &                                , ' land ( ', inbi0(ji),' x ', inbj0(ji),   & 
     866               & ' ), nb_points ', iszi0(ji)*iszj0(ji),' ( ', iszi0(ji),' x ', iszj0(ji),' )' 
     867         END DO 
     868         DEALLOCATE( inbi0, inbj0, iszi0, iszj0 ) 
     869         RETURN 
     870      ENDIF 
     871       
     872      DEALLOCATE( iszi0, iszj0 ) 
     873      inbij = inbijmax + 1        ! default: larger than possible 
     874      ii = isz0+1                 ! start from the end of the list (smaller subdomains) 
     875      DO WHILE( inbij > knbij )   ! while the number of ocean subdomains exceed the number of procs 
     876         ii = ii -1  
     877         ALLOCATE( llisoce(inbi0(ii), inbj0(ii)) ) 
     878         CALL mpp_init_isoce( inbi0(ii), inbj0(ii), llisoce )            ! must be done by all core 
     879         inbij = COUNT(llisoce) 
     880         DEALLOCATE( llisoce ) 
     881      END DO 
     882      knbi = inbi0(ii) 
     883      knbj = inbj0(ii) 
     884      DEALLOCATE( inbi0, inbj0 ) 
     885      ! 
     886   END SUBROUTINE mpp_init_bestpartition 
     887    
     888    
     889   SUBROUTINE mpp_init_landprop( propland ) 
     890      !!---------------------------------------------------------------------- 
     891      !!                  ***  ROUTINE mpp_init_landprop  *** 
     892      !! 
     893      !! ** Purpose : the the proportion of land points in the surface land-sea mask 
     894      !! 
     895      !! ** Method  : read iproc strips (of length jpiglo) of the land-sea mask 
     896      !!---------------------------------------------------------------------- 
     897      REAL(wp), INTENT(  out) :: propland    ! proportion of land points (between 0 and 1) 
     898      ! 
     899      INTEGER, DIMENSION(jpni*jpnj) ::   kusedom_1d 
     900      INTEGER :: inboce  
     901      INTEGER :: iproc, idiv, ijsz 
     902      INTEGER :: ijstr, ijend, ijcnt 
     903      LOGICAL, ALLOCATABLE, DIMENSION(:,:) ::   lloce 
     904      !!---------------------------------------------------------------------- 
     905      ! do nothing if there is no land-sea mask 
     906      IF( numbot == -1 .and. numbdy == -1 ) THEN 
     907         propland = 0. 
     908         RETURN 
     909      ENDIF 
     910 
     911      ! number of processes reading the bathymetry file  
     912      iproc = MINVAL( (/mppsize, jpjglo/2, 100/) )  ! read a least 2 lines, no more that 100 processes reading at the same time 
     913       
     914      ! we want to read iproc strips of the land-sea mask. -> pick up iproc processes among mppsize processes 
     915      IF( iproc == 1 ) THEN   ;   idiv = mppsize 
     916      ELSE                    ;   idiv = ( mppsize - 1 ) / ( iproc - 1 ) 
     917      ENDIF 
     918      ijsz = jpjglo / iproc 
     919      IF( narea <= MOD(jpjglo,iproc) ) ijsz = ijsz + 1 
     920       
     921      IF( MOD( narea-1, idiv ) == 0 .AND. (idiv /= 1 .OR. narea <= iproc ) ) THEN 
     922         ! 
     923         ijstr = (narea-1)*(jpjglo/iproc) + MIN(narea-1, MOD(jpjglo,iproc)) + 1 
     924         ijend = ijstr + ijsz - 1 
     925         ijcnt = ijend - ijstr + 1 
     926         ! 
     927         ALLOCATE( lloce(jpiglo, ijcnt) )   ! allocate the strip 
     928         CALL mpp_init_readbot_strip( ijstr, ijcnt, lloce ) 
     929         inboce = COUNT(lloce) 
     930         DEALLOCATE(lloce) 
     931         ! 
     932      ELSE 
     933         inboce = 0 
     934      ENDIF 
     935      CALL mpp_sum( 'mppini', inboce ) 
     936      ! 
     937      propland = REAL( jpiglo*jpjglo - inboce, wp ) / REAL( jpiglo*jpjglo, wp )  
     938      ! 
     939   END SUBROUTINE mpp_init_landprop 
     940    
     941    
     942   SUBROUTINE mpp_init_isoce( knbi, knbj, ldisoce ) 
     943      !!---------------------------------------------------------------------- 
     944      !!                  ***  ROUTINE mpp_init_nboce  *** 
     945      !! 
     946      !! ** Purpose : check for a mpi domain decomposition knbi x knbj which 
     947      !!              subdomains contain at least 1 ocean point 
     948      !! 
     949      !! ** Method  : read knbj strips (of length jpiglo) of the land-sea mask 
     950      !!---------------------------------------------------------------------- 
     951      INTEGER,                       INTENT(in   ) ::   knbi, knbj 
     952      LOGICAL, DIMENSION(knbi,knbj), INTENT(  out) ::   ldisoce   !  
     953      ! 
     954      INTEGER, DIMENSION(knbi,knbj) ::   inboce 
     955      INTEGER, DIMENSION(knbi*knbj) ::   inboce_1d 
     956      INTEGER :: idiv, i2read, inj 
     957      INTEGER :: iimax, ijmax 
     958      INTEGER :: ji,jj 
     959      LOGICAL, ALLOCATABLE, DIMENSION(:,:) ::   lloce 
     960      INTEGER, ALLOCATABLE, DIMENSION(:,:) ::   iimppt, ilci 
     961      INTEGER, ALLOCATABLE, DIMENSION(:,:) ::   ijmppt, ilcj 
     962      !!---------------------------------------------------------------------- 
     963      ! do nothing if there is no land-sea mask 
     964      IF( numbot == -1 .AND. numbdy == -1 ) THEN 
     965         ldisoce(:,:) = .TRUE. 
     966         RETURN 
     967      ENDIF 
     968 
     969      ! we want to read knbj strips of the land-sea mask. -> pick up knbj processes among mppsize processes 
     970      IF( knbj == 1 ) THEN   ;   idiv = mppsize 
     971      ELSE                   ;   idiv = ( mppsize - 1 ) / ( knbj - 1 ) 
     972      ENDIF 
     973      inboce(:,:) = 0 
     974      IF( MOD( narea-1, idiv ) == 0 .AND. (idiv /= 1 .OR. narea <= knbj ) ) THEN 
     975         ! 
     976         ALLOCATE( iimppt(knbi,knbj), ijmppt(knbi,knbj), ilci(knbi,knbj), ilcj(knbi,knbj) ) 
     977         CALL mpp_basic_decomposition( knbi, knbj, iimax, ijmax, iimppt, ijmppt, ilci, ilcj ) 
     978         ! 
     979         i2read = knbj / mppsize    ! strip number to be read by this process 
     980         IF( ( narea - 1 ) / idiv < MOD(knbj,mppsize) ) i2read = i2read + 1 
     981         DO jj = 1, i2read 
     982            ! strip number to be read (from 1 to knbj) 
     983            inj = ( narea - 1 ) * ( knbj / mppsize ) + MIN( MOD(knbj,mppsize), ( narea - 1 ) / idiv ) + jj 
     984            ALLOCATE( lloce(jpiglo, ilcj(1,inj)) )                              ! allocate the strip 
     985            CALL mpp_init_readbot_strip( ijmppt(1,inj), ilcj(1,inj), lloce )   ! read the strip 
     986            DO  ji = 1, knbi 
     987               inboce(ji,inj) = COUNT( lloce(iimppt(ji,1):iimppt(ji,1)+ilci(ji,1)-1,:) ) 
     988            END DO 
     989            DEALLOCATE(lloce) 
     990         END DO 
     991         ! 
     992         DEALLOCATE(iimppt, ijmppt, ilci, ilcj) 
     993         ! 
     994      ENDIF 
     995       
     996      inboce_1d = RESHAPE(inboce, (/ knbi*knbj /)) 
     997      CALL mpp_sum( 'mppini', inboce_1d ) 
     998      inboce = RESHAPE(inboce_1d, (/knbi, knbj/)) 
     999      ldisoce = inboce /= 0 
     1000      ! 
     1001   END SUBROUTINE mpp_init_isoce 
     1002    
     1003    
     1004   SUBROUTINE mpp_init_readbot_strip( kjstr, kjcnt, ldoce ) 
     1005      !!---------------------------------------------------------------------- 
     1006      !!                  ***  ROUTINE mpp_init_readbot_strip  *** 
     1007      !! 
     1008      !! ** Purpose : Read relevant bathymetric information in order to 
     1009      !!              provide a land/sea mask used for the elimination 
    6231010      !!              of land domains, in an mpp computation. 
    6241011      !! 
    625       !! ** Method  : Read the namelist ln_zco and ln_isfcav in namelist namzgr 
    626       !!              in order to choose the correct bathymetric information 
    627       !!              (file and variables)   
    628       !!---------------------------------------------------------------------- 
    629       INTEGER, DIMENSION(jpiglo,jpjglo), INTENT(out) ::   kmask   ! global domain  
    630    
    631       INTEGER :: inum   !: logical unit for configuration file 
    632       INTEGER :: ios    !: iostat error flag 
    633       INTEGER ::  ijstartrow                   ! temporary integers 
    634       REAL(wp), DIMENSION(jpiglo,jpjglo) ::   zbot, zbdy          ! global workspace 
    635       REAL(wp) ::   zidom , zjdom          ! local scalars 
    636       NAMELIST/nambdy/ ln_bdy, nb_bdy, ln_coords_file, cn_coords_file,           & 
    637            &             ln_mask_file, cn_mask_file, cn_dyn2d, nn_dyn2d_dta,     & 
    638            &             cn_dyn3d, nn_dyn3d_dta, cn_tra, nn_tra_dta,             &   
    639            &             ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp, rn_time_dmp_out, & 
    640            &             cn_ice, nn_ice_dta,                                     & 
    641            &             rn_ice_tem, rn_ice_sal, rn_ice_age,                     & 
    642            &             ln_vol, nn_volctl, nn_rimwidth, nb_jpk_bdy 
    643       !!---------------------------------------------------------------------- 
    644       ! 0. initialisation 
    645       ! ----------------- 
    646       CALL iom_open( cn_domcfg, inum ) 
    647       ! 
    648       ! ocean bottom level 
    649       CALL iom_get( inum, jpdom_unknown, 'bottom_level' , zbot , lrowattr=ln_use_jattr )  ! nb of ocean T-points 
    650       ! 
    651       CALL iom_close( inum ) 
    652       ! 
    653       ! 2D ocean mask (=1 if at least one level of the water column is ocean, =0 otherwise) 
    654       WHERE( zbot(:,:) > 0 )   ;   kmask(:,:) = 1 
    655       ELSEWHERE                ;   kmask(:,:) = 0 
    656       END WHERE 
    657    
    658       ! Adjust kmask with bdy_msk if it exists 
    659    
    660       REWIND( numnam_ref )              ! Namelist nambdy in reference namelist : BDY 
    661       READ  ( numnam_ref, nambdy, IOSTAT = ios, ERR = 903) 
    662 903   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy in reference namelist (mppini)', lwp ) 
    663       ! 
    664       REWIND( numnam_cfg )              ! Namelist nambdy in configuration namelist : BDY 
    665       READ  ( numnam_cfg, nambdy, IOSTAT = ios, ERR = 904 ) 
    666 904   IF( ios >  0 )   CALL ctl_nam ( ios , 'nambdy in configuration namelist (mppini)', lwp ) 
    667  
    668       IF( ln_bdy .AND. ln_mask_file ) THEN 
    669          CALL iom_open( cn_mask_file, inum ) 
    670          CALL iom_get ( inum, jpdom_unknown, 'bdy_msk', zbdy ) 
    671          CALL iom_close( inum ) 
    672          WHERE ( zbdy(:,:) <= 0. ) kmask = 0 
    673       ENDIF 
    674       ! 
    675    END SUBROUTINE mpp_init_mask 
     1012      !! ** Method  : read stipe of size (jpiglo,...) 
     1013      !!---------------------------------------------------------------------- 
     1014      INTEGER                         , INTENT(in   ) :: kjstr       !  
     1015      INTEGER                         , INTENT(in   ) :: kjcnt       !  
     1016      LOGICAL, DIMENSION(jpiglo,kjcnt), INTENT(  out) :: ldoce       !  
     1017      ! 
     1018      INTEGER                           ::   inumsave                     ! local logical unit 
     1019      REAL(wp), DIMENSION(jpiglo,kjcnt) ::   zbot, zbdy  
     1020      !!---------------------------------------------------------------------- 
     1021      ! 
     1022      inumsave = numout   ;   numout = numnul   !   redirect all print to /dev/null 
     1023      ! 
     1024      IF( numbot /= -1 ) THEN 
     1025         CALL iom_get( numbot, jpdom_unknown, 'bottom_level', zbot, kstart = (/1,kjstr/), kcount = (/jpiglo, kjcnt/) ) 
     1026      ELSE 
     1027         zbot(:,:) = 1.   ! put a non-null value 
     1028      ENDIF 
     1029 
     1030       IF( numbdy /= -1 ) THEN   ! Adjust  with bdy_msk if it exists     
     1031         CALL iom_get ( numbdy, jpdom_unknown, 'bdy_msk', zbdy, kstart = (/1,kjstr/), kcount = (/jpiglo, kjcnt/) ) 
     1032         zbot(:,:) = zbot(:,:) * zbdy(:,:) 
     1033      ENDIF 
     1034      ! 
     1035      ldoce = zbot > 0. 
     1036      numout = inumsave 
     1037      ! 
     1038   END SUBROUTINE mpp_init_readbot_strip 
    6761039 
    6771040 
     
    7201083      ! 
    7211084   END SUBROUTINE mpp_init_ioipsl   
    722  
    723  
    724    SUBROUTINE mpp_init_partition( num_pes ) 
    725       !!---------------------------------------------------------------------- 
    726       !!                 ***  ROUTINE mpp_init_partition  *** 
    727       !! 
    728       !! ** Purpose : 
    729       !! 
    730       !! ** Method  : 
    731       !!---------------------------------------------------------------------- 
    732       INTEGER, INTENT(in) ::   num_pes   ! The number of MPI processes we have 
    733       ! 
    734       INTEGER, PARAMETER :: nfactmax = 20 
    735       INTEGER :: nfact ! The no. of factors returned 
    736       INTEGER :: ierr  ! Error flag 
    737       INTEGER :: ji 
    738       INTEGER :: idiff, mindiff, imin ! For choosing pair of factors that are closest in value 
    739       INTEGER, DIMENSION(nfactmax) :: ifact ! Array of factors 
    740       !!---------------------------------------------------------------------- 
    741       ! 
    742       ierr = 0 
    743       ! 
    744       CALL factorise( ifact, nfactmax, nfact, num_pes, ierr ) 
    745       ! 
    746       IF( nfact <= 1 ) THEN 
    747          WRITE (numout, *) 'WARNING: factorisation of number of PEs failed' 
    748          WRITE (numout, *) '       : using grid of ',num_pes,' x 1' 
    749          jpnj = 1 
    750          jpni = num_pes 
    751       ELSE 
    752          ! Search through factors for the pair that are closest in value 
    753          mindiff = 1000000 
    754          imin    = 1 
    755          DO ji = 1, nfact-1, 2 
    756             idiff = ABS( ifact(ji) - ifact(ji+1) ) 
    757             IF( idiff < mindiff ) THEN 
    758                mindiff = idiff 
    759                imin = ji 
    760             ENDIF 
    761          END DO 
    762          jpnj = ifact(imin) 
    763          jpni = ifact(imin + 1) 
    764       ENDIF 
    765       ! 
    766       jpnij = jpni*jpnj 
    767       ! 
    768    END SUBROUTINE mpp_init_partition 
    769  
    770  
    771    SUBROUTINE factorise( kfax, kmaxfax, knfax, kn, kerr ) 
    772       !!---------------------------------------------------------------------- 
    773       !!                     ***  ROUTINE factorise  *** 
    774       !! 
    775       !! ** Purpose :   return the prime factors of n. 
    776       !!                knfax factors are returned in array kfax which is of 
    777       !!                maximum dimension kmaxfax. 
    778       !! ** Method  : 
    779       !!---------------------------------------------------------------------- 
    780       INTEGER                    , INTENT(in   ) ::   kn, kmaxfax 
    781       INTEGER                    , INTENT(  out) ::   kerr, knfax 
    782       INTEGER, DIMENSION(kmaxfax), INTENT(  out) ::   kfax 
    783       ! 
    784       INTEGER :: ifac, jl, inu 
    785       INTEGER, PARAMETER :: ntest = 14 
    786       INTEGER, DIMENSION(ntest) ::   ilfax 
    787       !!---------------------------------------------------------------------- 
    788       ! 
    789       ! lfax contains the set of allowed factors. 
    790       ilfax(:) = (/(2**jl,jl=ntest,1,-1)/) 
    791       ! 
    792       ! Clear the error flag and initialise output vars 
    793       kerr  = 0 
    794       kfax  = 1 
    795       knfax = 0 
    796       ! 
    797       IF( kn /= 1 ) THEN      ! Find the factors of n 
    798          ! 
    799          ! nu holds the unfactorised part of the number. 
    800          ! knfax holds the number of factors found. 
    801          ! l points to the allowed factor list. 
    802          ! ifac holds the current factor. 
    803          ! 
    804          inu   = kn 
    805          knfax = 0 
    806          ! 
    807          DO jl = ntest, 1, -1 
    808             ! 
    809             ifac = ilfax(jl) 
    810             IF( ifac > inu )   CYCLE 
    811             ! 
    812             ! Test whether the factor will divide. 
    813             ! 
    814             IF( MOD(inu,ifac) == 0 ) THEN 
    815                ! 
    816                knfax = knfax + 1            ! Add the factor to the list 
    817                IF( knfax > kmaxfax ) THEN 
    818                   kerr = 6 
    819                   write (*,*) 'FACTOR: insufficient space in factor array ', knfax 
    820                   return 
    821                ENDIF 
    822                kfax(knfax) = ifac 
    823                ! Store the other factor that goes with this one 
    824                knfax = knfax + 1 
    825                kfax(knfax) = inu / ifac 
    826                !WRITE (*,*) 'ARPDBG, factors ',knfax-1,' & ',knfax,' are ', kfax(knfax-1),' and ',kfax(knfax) 
    827             ENDIF 
    828             ! 
    829          END DO 
    830          ! 
    831       ENDIF 
    832       ! 
    833    END SUBROUTINE factorise 
    8341085 
    8351086 
     
    8961147   END SUBROUTINE mpp_init_nfdcom 
    8971148 
    898     
     1149 
    8991150#endif 
    9001151 
  • NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/nemogcm.F90

    r10297 r10330  
    234234#else 
    235235      IF    ( lk_oasis ) THEN   ;   CALL cpl_finalize   ! end coupling and mpp communications with OASIS 
    236       ELSEIF( lk_mpp   ) THEN   ;   CALL mppstop        ! end mpp communications 
     236      ELSEIF( lk_mpp   ) THEN   ;   CALL mppstop( ldfinal = .TRUE. )   ! end mpp communications 
    237237      ENDIF 
    238238#endif 
     
    349349         WRITE(numout,*) 
    350350         DO ji = 1, SIZE(cltxt) 
    351             IF( TRIM(cltxt (ji)) /= '' )   WRITE(numout,*) cltxt(ji)    ! control print of mynode 
     351            IF( TRIM(cltxt (ji)) /= '' )   WRITE(numout,*) TRIM(cltxt(ji))    ! control print of mynode 
    352352         END DO 
    353353         WRITE(numout,*) 
    354354         WRITE(numout,*) 
    355355         DO ji = 1, SIZE(cltxt2) 
    356             IF( TRIM(cltxt2(ji)) /= '' )   WRITE(numout,*) cltxt2(ji)   ! control print of domain size 
     356            IF( TRIM(cltxt2(ji)) /= '' )   WRITE(numout,*) TRIM(cltxt2(ji))   ! control print of domain size 
    357357         END DO 
    358358         ! 
     
    360360         ! 
    361361      ENDIF 
     362      ! open /dev/null file to be able to supress output write easily 
     363      CALL ctl_opn( numnul, '/dev/null', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE. ) 
     364      ! 
    362365      !                                      ! Domain decomposition 
    363366      CALL mpp_init                          ! MPP 
Note: See TracChangeset for help on using the changeset viewer.