New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 10425 for NEMO/trunk/src/OCE/LBC/mppini.F90 – NEMO

Ignore:
Timestamp:
2018-12-19T22:54:16+01:00 (5 years ago)
Author:
smasson
Message:

trunk: merge back dev_r10164_HPC09_ESIWACE_PREP_MERGE@10424 into the trunk

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk/src/OCE/LBC/mppini.F90

    r10068 r10425  
    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 ::   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                  !  
     151      INTEGER ::   inbi, inbj, iimax,  ijmax, icnt1, icnt2 
     152      LOGICAL ::   llbest 
    147153      INTEGER, ALLOCATABLE, DIMENSION(:)     ::   iin, ii_nono, ii_noea          ! 1D workspace 
    148154      INTEGER, ALLOCATABLE, DIMENSION(:)     ::   ijn, ii_noso, ii_nowe          !  -     - 
     
    151157      INTEGER, ALLOCATABLE, DIMENSION(:,:) ::   ilei, ildi, iono, ioea         !  -     - 
    152158      INTEGER, ALLOCATABLE, DIMENSION(:,:) ::   ilej, ildj, ioso, iowe         !  -     - 
    153       INTEGER, DIMENSION(jpiglo,jpjglo) ::   imask   ! 2D global domain workspace 
    154       !!---------------------------------------------------------------------- 
    155  
     159      LOGICAL, ALLOCATABLE, DIMENSION(:,:) ::   llisoce                        !  -     - 
     160      NAMELIST/nambdy/ ln_bdy, nb_bdy, ln_coords_file, cn_coords_file,           & 
     161           &             ln_mask_file, cn_mask_file, cn_dyn2d, nn_dyn2d_dta,     & 
     162           &             cn_dyn3d, nn_dyn3d_dta, cn_tra, nn_tra_dta,             &   
     163           &             ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp, rn_time_dmp_out, & 
     164           &             cn_ice, nn_ice_dta,                                     & 
     165           &             rn_ice_tem, rn_ice_sal, rn_ice_age,                     & 
     166           &             ln_vol, nn_volctl, nn_rimwidth, nb_jpk_bdy 
     167      !!---------------------------------------------------------------------- 
     168 
     169      ! do we need to take into account bdy_msk? 
     170      REWIND( numnam_ref )              ! Namelist nambdy in reference namelist : BDY 
     171      READ  ( numnam_ref, nambdy, IOSTAT = ios, ERR = 903) 
     172903   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy in reference namelist (mppini)', lwp ) 
     173      REWIND( numnam_cfg )              ! Namelist nambdy in configuration namelist : BDY 
     174      READ  ( numnam_cfg, nambdy, IOSTAT = ios, ERR = 904 ) 
     175904   IF( ios >  0 )   CALL ctl_nam ( ios , 'nambdy in configuration namelist (mppini)', lwp ) 
     176      ! 
     177      IF(               ln_read_cfg ) CALL iom_open( cn_domcfg,    numbot ) 
     178      IF( ln_bdy .AND. ln_mask_file ) CALL iom_open( cn_mask_file, numbdy ) 
     179      ! 
     180      !  1. Dimension arrays for subdomains 
     181      ! ----------------------------------- 
     182      ! 
    156183      ! If dimensions of processor grid weren't specified in the namelist file 
    157184      ! 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       ! 
     185      IF( jpni < 1 .OR. jpnj < 1 ) THEN 
     186         CALL mpp_init_bestpartition( mppsize, jpni, jpnj ) 
     187         llbest = .TRUE. 
     188      ELSE 
     189         CALL mpp_init_bestpartition( mppsize, inbi, inbj, icnt2 ) 
     190         CALL mpp_basic_decomposition( jpni, jpnj, jpimax, jpjmax ) 
     191         CALL mpp_basic_decomposition( inbi, inbj,  iimax,  ijmax ) 
     192         IF( iimax*ijmax < jpimax*jpjmax ) THEN 
     193            llbest = .FALSE. 
     194            icnt1 = jpni*jpnj - mppsize 
     195            WRITE(ctmp1,9000) '   The chosen domain decomposition ', jpni, ' x ', jpnj, ' with ', icnt1, ' land sub-domains' 
     196            WRITE(ctmp2,9000) '   has larger MPI subdomains (jpi = ', jpimax, ', jpj = ', jpjmax, ', jpi*jpj = ', jpimax*jpjmax, ')' 
     197            WRITE(ctmp3,9000) '   than the following domain decompostion ', inbi, ' x ', inbj, ' with ', icnt2, ' land sub-domains' 
     198            WRITE(ctmp4,9000) '   which MPI subdomains size is jpi = ', iimax, ', jpj = ', ijmax, ', jpi*jpj = ', iimax*ijmax, ' ' 
     199            CALL ctl_warn( 'mpp_init:', '~~~~~~~~ ', ctmp1, ctmp2, ctmp3, ctmp4, ' ', '    --- YOU ARE WASTING CPU... ---', ' ' ) 
     200         ELSE 
     201            llbest = .TRUE. 
     202         ENDIF 
     203      ENDIF 
     204       
     205      ! look for land mpi subdomains... 
     206      ALLOCATE( llisoce(jpni,jpnj) ) 
     207      CALL mpp_init_isoce( jpni, jpnj, llisoce ) 
     208      inijmin = COUNT( llisoce )   ! number of oce subdomains 
     209 
     210      IF( mppsize < inijmin ) THEN 
     211         WRITE(ctmp1,9001) '   With this specified domain decomposition: jpni = ', jpni, ' jpnj = ', jpnj 
     212         WRITE(ctmp2,9002) '   we can eliminate only ', jpni*jpnj - inijmin, ' land mpi subdomains therefore ' 
     213         WRITE(ctmp3,9001) '   the number of ocean mpi subdomains (', inijmin,') exceed the number of MPI processes:', mppsize 
     214         WRITE(ctmp4,*) '   ==>>> There is the list of best domain decompositions you should use: ' 
     215         CALL ctl_stop( 'mpp_init:', '~~~~~~~~ ', ctmp1, ctmp2, ctmp3, ctmp4 ) 
     216         CALL mpp_init_bestpartition( mppsize, ldlist = .TRUE. )   ! must be done by all core 
     217         CALL ctl_stop( 'STOP' ) 
     218      ENDIF 
     219 
     220      IF( mppsize > jpni*jpnj ) THEN 
     221         WRITE(ctmp1,9003) '   The number of mpi processes: ', mppsize 
     222         WRITE(ctmp2,9003) '   exceeds the maximum number of subdomains (ocean+land) = ', jpni*jpnj 
     223         WRITE(ctmp3,9001) '   defined by the following domain decomposition: jpni = ', jpni, ' jpnj = ', jpnj 
     224         WRITE(ctmp4,*) '   ==>>> There is the list of best domain decompositions you should use: ' 
     225         CALL ctl_stop( 'mpp_init:', '~~~~~~~~ ', ctmp1, ctmp2, ctmp3, ctmp4 ) 
     226         CALL mpp_init_bestpartition( mppsize, ldlist = .TRUE. )   ! must be done by all core 
     227         CALL ctl_stop( 'STOP' ) 
     228      ENDIF 
     229 
     230      jpnij = mppsize   ! force jpnij definition <-- remove as much land subdomains as needed to reach this condition 
     231      IF( mppsize > inijmin ) THEN 
     232         WRITE(ctmp1,9003) '   The number of mpi processes: ', mppsize 
     233         WRITE(ctmp2,9003) '   exceeds the maximum number of ocean subdomains = ', inijmin 
     234         WRITE(ctmp3,9002) '   we suppressed ', jpni*jpnj - mppsize, ' land subdomains ' 
     235         WRITE(ctmp4,9002) '   BUT we had to keep ', mppsize - inijmin, ' land subdomains that are useless...' 
     236         CALL ctl_warn( 'mpp_init:', '~~~~~~~~ ', ctmp1, ctmp2, ctmp3, ctmp4, ' ', '    --- YOU ARE WASTING CPU... ---', ' ' ) 
     237      ELSE   ! mppsize = inijmin 
     238         IF(lwp) THEN 
     239            IF(llbest) WRITE(numout,*) 'mpp_init: You use an optimal domain decomposition' 
     240            WRITE(numout,*) '~~~~~~~~ ' 
     241            WRITE(numout,9003) '   Number of mpi processes: ', mppsize 
     242            WRITE(numout,9003) '   Number of ocean subdomains = ', inijmin 
     243            WRITE(numout,9003) '   Number of suppressed land subdomains = ', jpni*jpnj - inijmin 
     244            WRITE(numout,*) 
     245         ENDIF 
     246      ENDIF 
     2479000  FORMAT (a, i4, a, i4, a, i7, a) 
     2489001  FORMAT (a, i4, a, i4) 
     2499002  FORMAT (a, i4, a) 
     2509003  FORMAT (a, i5) 
     251 
     252      IF( numbot /= -1 )   CALL iom_close( numbot ) 
     253      IF( numbdy /= -1 )   CALL iom_close( numbdy ) 
     254     
    164255      ALLOCATE(  nfiimpp(jpni,jpnj), nfipproc(jpni,jpnj), nfilcit(jpni,jpnj) ,    & 
    165256         &       nimppt(jpnij) , ibonit(jpnij) , nlcit(jpnij) , nlcjt(jpnij) ,    & 
     
    173264         &       ilej(jpni,jpnj), ildj(jpni,jpnj), ioso(jpni,jpnj), iowe(jpni,jpnj),   & 
    174265         &       STAT=ierr ) 
    175       CALL mpp_sum( ierr ) 
     266      CALL mpp_sum( 'mppini', ierr ) 
    176267      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'mpp_init: unable to allocate standard ocean arrays' ) 
    177268       
    178       ! 
    179269#if defined key_agrif 
    180270      IF( .NOT. Agrif_Root() ) THEN       ! AGRIF children: specific setting (cf. agrif_user.F90) 
     
    186276      ENDIF 
    187277#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 
     278      ! 
     279      !  2. Index arrays for subdomains 
    207280      ! ----------------------------------- 
    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. 
    212281      ! 
    213282      nreci = 2 * nn_hls 
    214283      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 ) 
     284      CALL mpp_basic_decomposition( jpni, jpnj, jpimax, jpjmax, iimppt, ijmppt, ilci, ilcj ) 
     285      nfiimpp(:,:) = iimppt(:,:) 
     286      nfilcit(:,:) = ilci(:,:) 
    236287      ! 
    237288      IF(lwp) THEN 
    238289         WRITE(numout,*) 
    239          WRITE(numout,*) 'mpp_init : MPI Message Passing MPI - domain lay out over processors' 
    240          WRITE(numout,*) '~~~~~~~~ ' 
     290         WRITE(numout,*) 'MPI Message Passing MPI - domain lay out over processors' 
     291         WRITE(numout,*) 
    241292         WRITE(numout,*) '   defines mpp subdomains' 
    242          WRITE(numout,*) '      iresti = ', iresti, ' jpni = ', jpni   
    243          WRITE(numout,*) '      irestj = ', irestj, ' jpnj = ', jpnj 
     293         WRITE(numout,*) '      jpni = ', jpni   
     294         WRITE(numout,*) '      jpnj = ', jpnj 
    244295         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  
     296         WRITE(numout,*) '      sum ilci(i,1) = ', sum(ilci(:,1)), ' jpiglo = ', jpiglo 
     297         WRITE(numout,*) '      sum ilcj(1,j) = ', sum(ilcj(1,:)), ' jpjglo = ', jpjglo 
     298      ENDIF 
     299      
    272300      ! 3. Subdomain description in the Regular Case 
    273301      ! -------------------------------------------- 
     
    277305      l_Jperio = jpnj == 1 .AND. (jperio == 2 .OR. jperio == 7) 
    278306       
    279       icont = -1 
    280307      DO jarea = 1, jpni*jpnj 
     308         ! 
    281309         iarea0 = jarea - 1 
    282310         ii = 1 + MOD(iarea0,jpni) 
     
    334362         ENDIF 
    335363         ! 
    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 
     364      END DO 
     365 
     366      ! 4. deal with land subdomains 
     367      ! ---------------------------- 
     368      ! 
     369      ! specify which subdomains are oce subdomains; other are land subdomains 
     370      ipproc(:,:) = -1 
     371      icont = -1 
     372      DO jarea = 1, jpni*jpnj 
     373         iarea0 = jarea - 1 
     374         ii = 1 + MOD(iarea0,jpni) 
     375         ij = 1 +     iarea0/jpni 
     376         IF( llisoce(ii,ij) ) THEN 
    345377            icont = icont + 1 
    346378            ipproc(ii,ij) = icont 
     
    349381         ENDIF 
    350382      END DO 
    351       ! 
     383      ! if needed add some land subdomains to reach jpnij active subdomains 
     384      i2add = jpnij - inijmin 
     385      DO jarea = 1, jpni*jpnj 
     386         iarea0 = jarea - 1 
     387         ii = 1 + MOD(iarea0,jpni) 
     388         ij = 1 +     iarea0/jpni 
     389         IF( .NOT. llisoce(ii,ij) .AND. i2add > 0 ) THEN 
     390            icont = icont + 1 
     391            ipproc(ii,ij) = icont 
     392            iin(icont+1) = ii 
     393            ijn(icont+1) = ij 
     394            i2add = i2add - 1 
     395         ENDIF 
     396      END DO 
    352397      nfipproc(:,:) = ipproc(:,:) 
    353398 
    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 
     399      ! neighbour treatment: change ibondi, ibondj if next to a land zone 
     400      DO jarea = 1, jpni*jpnj 
     401         ii = 1 + MOD( jarea-1  , jpni ) 
     402         ij = 1 +     (jarea-1) / jpni 
     403         ! land-only area with an active n neigbour 
     404         IF ( ipproc(ii,ij) == -1 .AND. 0 <= iono(ii,ij) .AND. iono(ii,ij) <= jpni*jpnj-1 ) THEN 
     405            iino = 1 + MOD( iono(ii,ij) , jpni )                    ! ii index of this n neigbour 
     406            ijno = 1 +      iono(ii,ij) / jpni                      ! ij index of this n neigbour 
     407            ! In case of north fold exchange: I am the n neigbour of my n neigbour!! (#1057) 
     408            ! --> for northern neighbours of northern row processors (in case of north-fold) 
     409            !     need to reverse the LOGICAL direction of communication  
     410            idir = 1                                           ! we are indeed the s neigbour of this n neigbour 
     411            IF( ij == jpnj .AND. ijno == jpnj )   idir = -1    ! both are on the last row, we are in fact the n neigbour 
     412            IF( ibondj(iino,ijno) == idir     )   ibondj(iino,ijno) =   2     ! this n neigbour had only a s/n neigbour -> no more 
     413            IF( ibondj(iino,ijno) == 0        )   ibondj(iino,ijno) = -idir   ! this n neigbour had both, n-s neighbours -> keep 1 
     414         ENDIF 
     415         ! land-only area with an active s neigbour 
     416         IF( ipproc(ii,ij) == -1 .AND. 0 <= ioso(ii,ij) .AND. ioso(ii,ij) <= jpni*jpnj-1 ) THEN 
     417            iiso = 1 + MOD( ioso(ii,ij) , jpni )                    ! ii index of this s neigbour 
     418            ijso = 1 +      ioso(ii,ij) / jpni                      ! ij index of this s neigbour 
     419            IF( ibondj(iiso,ijso) == -1 )   ibondj(iiso,ijso) = 2   ! this s neigbour had only a n neigbour    -> no more neigbour 
     420            IF( ibondj(iiso,ijso) ==  0 )   ibondj(iiso,ijso) = 1   ! this s neigbour had both, n-s neighbours -> keep s neigbour 
     421         ENDIF 
     422         ! land-only area with an active e neigbour 
     423         IF( ipproc(ii,ij) == -1 .AND. 0 <= ioea(ii,ij) .AND. ioea(ii,ij) <= jpni*jpnj-1 ) THEN 
     424            iiea = 1 + MOD( ioea(ii,ij) , jpni )                    ! ii index of this e neigbour 
     425            ijea = 1 +      ioea(ii,ij) / jpni                      ! ij index of this e neigbour 
     426            IF( ibondi(iiea,ijea) == 1 )   ibondi(iiea,ijea) =  2   ! this e neigbour had only a w neigbour    -> no more neigbour 
     427            IF( ibondi(iiea,ijea) == 0 )   ibondi(iiea,ijea) = -1   ! this e neigbour had both, e-w neighbours -> keep e neigbour 
     428         ENDIF 
     429         ! land-only area with an active w neigbour 
     430         IF( ipproc(ii,ij) == -1 .AND. 0 <= iowe(ii,ij) .AND. iowe(ii,ij) <= jpni*jpnj-1) THEN 
     431            iiwe = 1 + MOD( iowe(ii,ij) , jpni )                    ! ii index of this w neigbour 
     432            ijwe = 1 +      iowe(ii,ij) / jpni                      ! ij index of this w neigbour 
     433            IF( ibondi(iiwe,ijwe) == -1 )   ibondi(iiwe,ijwe) = 2   ! this w neigbour had only a e neigbour    -> no more neigbour 
     434            IF( ibondi(iiwe,ijwe) ==  0 )   ibondi(iiwe,ijwe) = 1   ! this w neigbour had both, e-w neighbours -> keep w neigbour 
     435         ENDIF 
     436      END DO 
     437 
     438      ! Update il[de][ij] according to modified ibond[ij] 
     439      ! ---------------------- 
     440      DO jproc = 1, jpnij 
     441         ii = iin(jproc) 
     442         ij = ijn(jproc) 
     443         IF( ibondi(ii,ij) == -1 .OR. ibondi(ii,ij) == 2 ) ildi(ii,ij) =  1 
     444         IF( ibondi(ii,ij) ==  1 .OR. ibondi(ii,ij) == 2 ) ilei(ii,ij) = ilci(ii,ij) 
     445         IF( ibondj(ii,ij) == -1 .OR. ibondj(ii,ij) == 2 ) ildj(ii,ij) =  1 
     446         IF( ibondj(ii,ij) ==  1 .OR. ibondj(ii,ij) == 2 ) ilej(ii,ij) = ilcj(ii,ij) 
     447      END DO 
     448       
     449      ! 5. Subdomain print 
    363450      ! ------------------ 
    364451      IF(lwp) THEN 
     
    385472 9404    FORMAT('           *  '   ,20('      ',i3,'   *   ') ) 
    386473      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 
    438474          
    439475      ! just to save nono etc for all proc 
     
    516552         njmppt(jproc) = ijmppt(ii,ij)  
    517553      END DO 
    518       nfilcit(:,:) = ilci(:,:) 
    519554 
    520555      ! Save processor layout in ascii file 
     
    610645         &       iimppt, ijmppt, ibondi, ibondj, ipproc, ipolj,   & 
    611646         &       ilci, ilcj, ilei, ilej, ildi, ildj,              & 
    612          &       iono, ioea, ioso, iowe) 
     647         &       iono, ioea, ioso, iowe, llisoce) 
    613648      ! 
    614649    END SUBROUTINE mpp_init 
    615650 
    616651 
    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 
     652    SUBROUTINE mpp_basic_decomposition( knbi, knbj, kimax, kjmax, kimppt, kjmppt, klci, klcj) 
     653      !!---------------------------------------------------------------------- 
     654      !!                  ***  ROUTINE mpp_basic_decomposition  *** 
     655      !!                     
     656      !! ** Purpose :   Lay out the global domain over processors. 
     657      !! 
     658      !! ** Method  :   Global domain is distributed in smaller local domains. 
     659      !! 
     660      !! ** Action : - set for all knbi*knbj domains: 
     661      !!                    kimppt     : longitudinal index 
     662      !!                    kjmppt     : latitudinal  index 
     663      !!                    klci       : first dimension 
     664      !!                    klcj       : second dimension 
     665      !!---------------------------------------------------------------------- 
     666      INTEGER,                                 INTENT(in   ) ::   knbi, knbj 
     667      INTEGER,                                 INTENT(  out) ::   kimax, kjmax 
     668      INTEGER, DIMENSION(knbi,knbj), OPTIONAL, INTENT(  out) ::   kimppt, kjmppt 
     669      INTEGER, DIMENSION(knbi,knbj), OPTIONAL, INTENT(  out) ::   klci, klcj 
     670      ! 
     671      INTEGER ::   ji, jj 
     672      INTEGER ::   iresti, irestj 
     673      INTEGER ::   ireci, irecj 
     674      !!---------------------------------------------------------------------- 
     675      ! 
     676#if defined key_nemocice_decomp 
     677      kimax = ( nx_global+2-2*nn_hls + (knbi-1) ) / knbi + 2*nn_hls    ! first  dim. 
     678      kjmax = ( ny_global+2-2*nn_hls + (knbj-1) ) / knbj + 2*nn_hls    ! second dim.  
     679#else 
     680      kimax = ( jpiglo - 2*nn_hls + (knbi-1) ) / knbi + 2*nn_hls    ! first  dim. 
     681      kjmax = ( jpjglo - 2*nn_hls + (knbj-1) ) / knbj + 2*nn_hls    ! second dim. 
     682#endif 
     683      IF( .NOT. PRESENT(kimppt) ) RETURN 
     684      ! 
     685      !  1. Dimension arrays for subdomains 
     686      ! ----------------------------------- 
     687      !  Computation of local domain sizes klci() klcj() 
     688      !  These dimensions depend on global sizes knbi,knbj and jpiglo,jpjglo 
     689      !  The subdomains are squares lesser than or equal to the global 
     690      !  dimensions divided by the number of processors minus the overlap array. 
     691      ! 
     692      ireci = 2 * nn_hls 
     693      irecj = 2 * nn_hls 
     694      iresti = 1 + MOD( jpiglo - ireci -1 , knbi ) 
     695      irestj = 1 + MOD( jpjglo - irecj -1 , knbj ) 
     696      ! 
     697      !  Need to use kimax and kjmax here since jpi and jpj not yet defined 
     698#if defined key_nemocice_decomp 
     699      ! Change padding to be consistent with CICE 
     700      klci(1:knbi-1      ,:) = kimax 
     701      klci(knbi          ,:) = jpiglo - (knbi - 1) * (kimax - nreci) 
     702      klcj(:,      1:knbj-1) = kjmax 
     703      klcj(:,          knbj) = jpjglo - (knbj - 1) * (kjmax - nrecj) 
     704#else 
     705      klci(1:iresti      ,:) = kimax 
     706      klci(iresti+1:knbi ,:) = kimax-1 
     707      klcj(:,      1:irestj) = kjmax 
     708      klcj(:, irestj+1:knbj) = kjmax-1 
     709#endif 
     710 
     711      !  2. Index arrays for subdomains 
     712      ! ------------------------------- 
     713      kimppt(:,:) = 1 
     714      kjmppt(:,:) = 1 
     715      ! 
     716      IF( knbi > 1 ) THEN 
     717         DO jj = 1, knbj 
     718            DO ji = 2, knbi 
     719               kimppt(ji,jj) = kimppt(ji-1,jj) + klci(ji-1,jj) - ireci 
     720            END DO 
     721         END DO 
     722      ENDIF 
     723      ! 
     724      IF( knbj > 1 )THEN 
     725         DO jj = 2, knbj 
     726            DO ji = 1, knbi 
     727               kjmppt(ji,jj) = kjmppt(ji,jj-1) + klcj(ji,jj-1) - irecj 
     728            END DO 
     729         END DO 
     730      ENDIF 
     731       
     732   END SUBROUTINE mpp_basic_decomposition 
     733 
     734 
     735   SUBROUTINE mpp_init_bestpartition( knbij, knbi, knbj, knbcnt, ldlist ) 
     736      !!---------------------------------------------------------------------- 
     737      !!                 ***  ROUTINE mpp_init_bestpartition  *** 
     738      !! 
     739      !! ** Purpose : 
     740      !! 
     741      !! ** Method  : 
     742      !!---------------------------------------------------------------------- 
     743      INTEGER,           INTENT(in   ) ::   knbij         ! total number if subdomains               (knbi*knbj) 
     744      INTEGER, OPTIONAL, INTENT(  out) ::   knbi, knbj    ! number if subdomains along i and j (knbi and knbj) 
     745      INTEGER, OPTIONAL, INTENT(  out) ::   knbcnt        ! number of land subdomains 
     746      LOGICAL, OPTIONAL, INTENT(in   ) ::   ldlist        ! .true.: print the list the best domain decompositions (with land) 
     747      ! 
     748      INTEGER :: ji, jj, ii, iitarget 
     749      INTEGER :: iszitst, iszjtst 
     750      INTEGER :: isziref, iszjref 
     751      INTEGER :: inbij, iszij 
     752      INTEGER :: inbimax, inbjmax, inbijmax 
     753      INTEGER :: isz0, isz1 
     754      INTEGER, DIMENSION(  :), ALLOCATABLE :: indexok 
     755      INTEGER, DIMENSION(  :), ALLOCATABLE :: inbi0, inbj0, inbij0   ! number of subdomains along i,j 
     756      INTEGER, DIMENSION(  :), ALLOCATABLE :: iszi0, iszj0, iszij0   ! max size of the subdomains along i,j 
     757      INTEGER, DIMENSION(  :), ALLOCATABLE :: inbi1, inbj1, inbij1   ! number of subdomains along i,j 
     758      INTEGER, DIMENSION(  :), ALLOCATABLE :: iszi1, iszj1, iszij1   ! max size of the subdomains along i,j 
     759      LOGICAL :: llist 
     760      LOGICAL, DIMENSION(:,:), ALLOCATABLE :: llmsk2d                 ! max size of the subdomains along i,j 
     761      LOGICAL, DIMENSION(:,:), ALLOCATABLE :: llisoce              !  -     - 
     762      REAL(wp)::   zpropland 
     763      !!---------------------------------------------------------------------- 
     764      ! 
     765      llist = .FALSE. 
     766      IF( PRESENT(ldlist) ) llist = ldlist 
     767 
     768      CALL mpp_init_landprop( zpropland )                      ! get the proportion of land point over the gloal domain 
     769      inbij = NINT( REAL(knbij, wp) / ( 1.0 - zpropland ) )    ! define the largest possible value for jpni*jpnj 
     770      ! 
     771      IF( llist ) THEN   ;   inbijmax = inbij*2 
     772      ELSE               ;   inbijmax = inbij 
     773      ENDIF 
     774      ! 
     775      ALLOCATE(inbi0(inbijmax),inbj0(inbijmax),iszi0(inbijmax),iszj0(inbijmax)) 
     776      ! 
     777      inbimax = 0 
     778      inbjmax = 0 
     779      isziref = jpiglo*jpjglo+1 
     780      iszjref = jpiglo*jpjglo+1 
     781      ! 
     782      ! get the list of knbi that gives a smaller jpimax than knbi-1 
     783      ! get the list of knbj that gives a smaller jpjmax than knbj-1 
     784      DO ji = 1, inbijmax       
     785#if defined key_nemocice_decomp 
     786         iszitst = ( nx_global+2-2*nn_hls + (ji-1) ) / ji + 2*nn_hls    ! first  dim. 
     787#else 
     788         iszitst = ( jpiglo - 2*nn_hls + (ji-1) ) / ji + 2*nn_hls 
     789#endif 
     790         IF( iszitst < isziref ) THEN 
     791            isziref = iszitst 
     792            inbimax = inbimax + 1 
     793            inbi0(inbimax) = ji 
     794            iszi0(inbimax) = isziref 
     795         ENDIF 
     796#if defined key_nemocice_decomp 
     797         iszjtst = ( ny_global+2-2*nn_hls + (ji-1) ) / ji + 2*nn_hls    ! first  dim. 
     798#else 
     799         iszjtst = ( jpjglo - 2*nn_hls + (ji-1) ) / ji + 2*nn_hls 
     800#endif 
     801         IF( iszjtst < iszjref ) THEN 
     802            iszjref = iszjtst 
     803            inbjmax = inbjmax + 1 
     804            inbj0(inbjmax) = ji 
     805            iszj0(inbjmax) = iszjref 
     806         ENDIF 
     807      END DO 
     808 
     809      ! combine these 2 lists to get all possible knbi*knbj <  inbijmax 
     810      ALLOCATE( llmsk2d(inbimax,inbjmax) ) 
     811      DO jj = 1, inbjmax 
     812         DO ji = 1, inbimax 
     813            IF ( inbi0(ji) * inbj0(jj) <= inbijmax ) THEN   ;   llmsk2d(ji,jj) = .TRUE. 
     814            ELSE                                            ;   llmsk2d(ji,jj) = .FALSE. 
     815            ENDIF 
     816         END DO 
     817      END DO 
     818      isz1 = COUNT(llmsk2d) 
     819      ALLOCATE( inbi1(isz1), inbj1(isz1), iszi1(isz1), iszj1(isz1) ) 
     820      ii = 0 
     821      DO jj = 1, inbjmax 
     822         DO ji = 1, inbimax 
     823            IF( llmsk2d(ji,jj) .EQV. .TRUE. ) THEN 
     824               ii = ii + 1 
     825               inbi1(ii) = inbi0(ji) 
     826               inbj1(ii) = inbj0(jj) 
     827               iszi1(ii) = iszi0(ji) 
     828               iszj1(ii) = iszj0(jj) 
     829            END IF 
     830         END DO 
     831      END DO 
     832      DEALLOCATE( inbi0, inbj0, iszi0, iszj0 ) 
     833      DEALLOCATE( llmsk2d ) 
     834 
     835      ALLOCATE( inbij1(isz1), iszij1(isz1) ) 
     836      inbij1(:) = inbi1(:) * inbj1(:) 
     837      iszij1(:) = iszi1(:) * iszj1(:) 
     838 
     839      ! if therr is no land and no print 
     840      IF( .NOT. llist .AND. numbot == -1 .AND. numbdy == -1 ) THEN 
     841         ! get the smaller partition which gives the smallest subdomain size 
     842         ii = MINLOC(inbij1, mask = iszij1 == MINVAL(iszij1), dim = 1) 
     843         knbi = inbi1(ii) 
     844         knbj = inbj1(ii) 
     845         IF(PRESENT(knbcnt))   knbcnt = 0 
     846         DEALLOCATE( inbi1, inbj1, inbij1, iszi1, iszj1, iszij1 ) 
     847         RETURN 
     848      ENDIF 
     849 
     850      ! extract only the partitions which reduce the subdomain size in comparison with smaller partitions 
     851      ALLOCATE( indexok(isz1) )                                 ! to store indices of the best partitions 
     852      isz0 = 0                                                  ! number of best partitions      
     853      inbij = 1                                                 ! start with the min value of inbij1 => 1 
     854      iszij = jpiglo*jpjglo+1                                   ! default: larger than global domain 
     855      DO WHILE( inbij <= inbijmax )                             ! if we did not reach the max of inbij1 
     856         ii = MINLOC(iszij1, mask = inbij1 == inbij, dim = 1)   ! warning: send back the first occurence if multiple results 
     857         IF ( iszij1(ii) < iszij ) THEN 
     858            isz0 = isz0 + 1 
     859            indexok(isz0) = ii 
     860            iszij = iszij1(ii) 
     861         ENDIF 
     862         inbij = MINVAL(inbij1, mask = inbij1 > inbij)   ! warning: return largest integer value if mask = .false. everywhere 
     863      END DO 
     864      DEALLOCATE( inbij1, iszij1 ) 
     865 
     866      ! keep only the best partitions (sorted by increasing order of subdomains number and decreassing subdomain size) 
     867      ALLOCATE( inbi0(isz0), inbj0(isz0), iszi0(isz0), iszj0(isz0) ) 
     868      DO ji = 1, isz0 
     869         ii = indexok(ji) 
     870         inbi0(ji) = inbi1(ii) 
     871         inbj0(ji) = inbj1(ii) 
     872         iszi0(ji) = iszi1(ii) 
     873         iszj0(ji) = iszj1(ii) 
     874      END DO 
     875      DEALLOCATE( indexok, inbi1, inbj1, iszi1, iszj1 ) 
     876 
     877      IF( llist ) THEN  ! we print about 21 best partitions 
     878         IF(lwp) THEN 
     879            WRITE(numout,*) 
     880            WRITE(numout,         *) '                  For your information:' 
     881            WRITE(numout,'(a,i5,a)') '  list of the best partitions around ',   knbij, ' mpi processes' 
     882            WRITE(numout,         *) '  --------------------------------------', '-----', '--------------' 
     883            WRITE(numout,*) 
     884         END IF 
     885         iitarget = MINLOC( inbi0(:)*inbj0(:), mask = inbi0(:)*inbj0(:) >= knbij, dim = 1 ) 
     886         DO ji = MAX(1,iitarget-10), MIN(isz0,iitarget+10) 
     887            ALLOCATE( llisoce(inbi0(ji), inbj0(ji)) ) 
     888            CALL mpp_init_isoce( inbi0(ji), inbj0(ji), llisoce ) ! Warning: must be call by all cores (call mpp_sum) 
     889            inbij = COUNT(llisoce) 
     890            DEALLOCATE( llisoce ) 
     891            IF(lwp) WRITE(numout,'(a, i5, a, i5, a, i4, a, i4, a, i9, a, i5, a, i5, a)')    & 
     892               &     'nb_cores ' , inbij,' oce + ', inbi0(ji)*inbj0(ji) - inbij             & 
     893               &                                , ' land ( ', inbi0(ji),' x ', inbj0(ji),   & 
     894               & ' ), nb_points ', iszi0(ji)*iszj0(ji),' ( ', iszi0(ji),' x ', iszj0(ji),' )' 
     895         END DO 
     896         DEALLOCATE( inbi0, inbj0, iszi0, iszj0 ) 
     897         RETURN 
     898      ENDIF 
     899       
     900      DEALLOCATE( iszi0, iszj0 ) 
     901      inbij = inbijmax + 1        ! default: larger than possible 
     902      ii = isz0+1                 ! start from the end of the list (smaller subdomains) 
     903      DO WHILE( inbij > knbij )   ! while the number of ocean subdomains exceed the number of procs 
     904         ii = ii -1  
     905         ALLOCATE( llisoce(inbi0(ii), inbj0(ii)) ) 
     906         CALL mpp_init_isoce( inbi0(ii), inbj0(ii), llisoce )            ! must be done by all core 
     907         inbij = COUNT(llisoce) 
     908         DEALLOCATE( llisoce ) 
     909      END DO 
     910      knbi = inbi0(ii) 
     911      knbj = inbj0(ii) 
     912      IF(PRESENT(knbcnt))   knbcnt = knbi * knbj - inbij 
     913      DEALLOCATE( inbi0, inbj0 ) 
     914      ! 
     915   END SUBROUTINE mpp_init_bestpartition 
     916    
     917    
     918   SUBROUTINE mpp_init_landprop( propland ) 
     919      !!---------------------------------------------------------------------- 
     920      !!                  ***  ROUTINE mpp_init_landprop  *** 
     921      !! 
     922      !! ** Purpose : the the proportion of land points in the surface land-sea mask 
     923      !! 
     924      !! ** Method  : read iproc strips (of length jpiglo) of the land-sea mask 
     925      !!---------------------------------------------------------------------- 
     926      REAL(wp), INTENT(  out) :: propland    ! proportion of land points (between 0 and 1) 
     927      ! 
     928      INTEGER, DIMENSION(jpni*jpnj) ::   kusedom_1d 
     929      INTEGER :: inboce  
     930      INTEGER :: iproc, idiv, ijsz 
     931      INTEGER :: ijstr, ijend, ijcnt 
     932      LOGICAL, ALLOCATABLE, DIMENSION(:,:) ::   lloce 
     933      !!---------------------------------------------------------------------- 
     934      ! do nothing if there is no land-sea mask 
     935      IF( numbot == -1 .and. numbdy == -1 ) THEN 
     936         propland = 0. 
     937         RETURN 
     938      ENDIF 
     939 
     940      ! number of processes reading the bathymetry file  
     941      iproc = MINVAL( (/mppsize, jpjglo/2, 100/) )  ! read a least 2 lines, no more that 100 processes reading at the same time 
     942       
     943      ! we want to read iproc strips of the land-sea mask. -> pick up iproc processes among mppsize processes 
     944      IF( iproc == 1 ) THEN   ;   idiv = mppsize 
     945      ELSE                    ;   idiv = ( mppsize - 1 ) / ( iproc - 1 ) 
     946      ENDIF 
     947      ijsz = jpjglo / iproc 
     948      IF( narea <= MOD(jpjglo,iproc) ) ijsz = ijsz + 1 
     949       
     950      IF( MOD( narea-1, idiv ) == 0 .AND. (idiv /= 1 .OR. narea <= iproc ) ) THEN 
     951         ! 
     952         ijstr = (narea-1)*(jpjglo/iproc) + MIN(narea-1, MOD(jpjglo,iproc)) + 1 
     953         ijend = ijstr + ijsz - 1 
     954         ijcnt = ijend - ijstr + 1 
     955         ! 
     956         ALLOCATE( lloce(jpiglo, ijcnt) )   ! allocate the strip 
     957         CALL mpp_init_readbot_strip( ijstr, ijcnt, lloce ) 
     958         inboce = COUNT(lloce) 
     959         DEALLOCATE(lloce) 
     960         ! 
     961      ELSE 
     962         inboce = 0 
     963      ENDIF 
     964      CALL mpp_sum( 'mppini', inboce ) 
     965      ! 
     966      propland = REAL( jpiglo*jpjglo - inboce, wp ) / REAL( jpiglo*jpjglo, wp )  
     967      ! 
     968   END SUBROUTINE mpp_init_landprop 
     969    
     970    
     971   SUBROUTINE mpp_init_isoce( knbi, knbj, ldisoce ) 
     972      !!---------------------------------------------------------------------- 
     973      !!                  ***  ROUTINE mpp_init_nboce  *** 
     974      !! 
     975      !! ** Purpose : check for a mpi domain decomposition knbi x knbj which 
     976      !!              subdomains contain at least 1 ocean point 
     977      !! 
     978      !! ** Method  : read knbj strips (of length jpiglo) of the land-sea mask 
     979      !!---------------------------------------------------------------------- 
     980      INTEGER,                       INTENT(in   ) ::   knbi, knbj 
     981      LOGICAL, DIMENSION(knbi,knbj), INTENT(  out) ::   ldisoce   !  
     982      ! 
     983      INTEGER, DIMENSION(knbi,knbj) ::   inboce 
     984      INTEGER, DIMENSION(knbi*knbj) ::   inboce_1d 
     985      INTEGER :: idiv, i2read, inj 
     986      INTEGER :: iimax, ijmax 
     987      INTEGER :: ji,jj 
     988      LOGICAL, ALLOCATABLE, DIMENSION(:,:) ::   lloce 
     989      INTEGER, ALLOCATABLE, DIMENSION(:,:) ::   iimppt, ilci 
     990      INTEGER, ALLOCATABLE, DIMENSION(:,:) ::   ijmppt, ilcj 
     991      !!---------------------------------------------------------------------- 
     992      ! do nothing if there is no land-sea mask 
     993      IF( numbot == -1 .AND. numbdy == -1 ) THEN 
     994         ldisoce(:,:) = .TRUE. 
     995         RETURN 
     996      ENDIF 
     997 
     998      ! we want to read knbj strips of the land-sea mask. -> pick up knbj processes among mppsize processes 
     999      IF( knbj == 1 ) THEN   ;   idiv = mppsize 
     1000      ELSE                   ;   idiv = ( mppsize - 1 ) / ( knbj - 1 ) 
     1001      ENDIF 
     1002      inboce(:,:) = 0 
     1003      IF( MOD( narea-1, idiv ) == 0 .AND. (idiv /= 1 .OR. narea <= knbj ) ) THEN 
     1004         ! 
     1005         ALLOCATE( iimppt(knbi,knbj), ijmppt(knbi,knbj), ilci(knbi,knbj), ilcj(knbi,knbj) ) 
     1006         CALL mpp_basic_decomposition( knbi, knbj, iimax, ijmax, iimppt, ijmppt, ilci, ilcj ) 
     1007         ! 
     1008         i2read = knbj / mppsize    ! strip number to be read by this process 
     1009         IF( ( narea - 1 ) / idiv < MOD(knbj,mppsize) ) i2read = i2read + 1 
     1010         DO jj = 1, i2read 
     1011            ! strip number to be read (from 1 to knbj) 
     1012            inj = ( narea - 1 ) * ( knbj / mppsize ) + MIN( MOD(knbj,mppsize), ( narea - 1 ) / idiv ) + jj 
     1013            ALLOCATE( lloce(jpiglo, ilcj(1,inj)) )                              ! allocate the strip 
     1014            CALL mpp_init_readbot_strip( ijmppt(1,inj), ilcj(1,inj), lloce )   ! read the strip 
     1015            DO  ji = 1, knbi 
     1016               inboce(ji,inj) = COUNT( lloce(iimppt(ji,1):iimppt(ji,1)+ilci(ji,1)-1,:) ) 
     1017            END DO 
     1018            DEALLOCATE(lloce) 
     1019         END DO 
     1020         ! 
     1021         DEALLOCATE(iimppt, ijmppt, ilci, ilcj) 
     1022         ! 
     1023      ENDIF 
     1024       
     1025      inboce_1d = RESHAPE(inboce, (/ knbi*knbj /)) 
     1026      CALL mpp_sum( 'mppini', inboce_1d ) 
     1027      inboce = RESHAPE(inboce_1d, (/knbi, knbj/)) 
     1028      ldisoce = inboce /= 0 
     1029      ! 
     1030   END SUBROUTINE mpp_init_isoce 
     1031    
     1032    
     1033   SUBROUTINE mpp_init_readbot_strip( kjstr, kjcnt, ldoce ) 
     1034      !!---------------------------------------------------------------------- 
     1035      !!                  ***  ROUTINE mpp_init_readbot_strip  *** 
     1036      !! 
     1037      !! ** Purpose : Read relevant bathymetric information in order to 
     1038      !!              provide a land/sea mask used for the elimination 
    6231039      !!              of land domains, in an mpp computation. 
    6241040      !! 
    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 
     1041      !! ** Method  : read stipe of size (jpiglo,...) 
     1042      !!---------------------------------------------------------------------- 
     1043      INTEGER                         , INTENT(in   ) :: kjstr       !  
     1044      INTEGER                         , INTENT(in   ) :: kjcnt       !  
     1045      LOGICAL, DIMENSION(jpiglo,kjcnt), INTENT(  out) :: ldoce       !  
     1046      ! 
     1047      INTEGER                           ::   inumsave                     ! local logical unit 
     1048      REAL(wp), DIMENSION(jpiglo,kjcnt) ::   zbot, zbdy  
     1049      !!---------------------------------------------------------------------- 
     1050      ! 
     1051      inumsave = numout   ;   numout = numnul   !   redirect all print to /dev/null 
     1052      ! 
     1053      IF( numbot /= -1 ) THEN 
     1054         CALL iom_get( numbot, jpdom_unknown, 'bottom_level', zbot, kstart = (/1,kjstr/), kcount = (/jpiglo, kjcnt/) ) 
     1055      ELSE 
     1056         zbot(:,:) = 1.   ! put a non-null value 
     1057      ENDIF 
     1058 
     1059       IF( numbdy /= -1 ) THEN   ! Adjust  with bdy_msk if it exists     
     1060         CALL iom_get ( numbdy, jpdom_unknown, 'bdy_msk', zbdy, kstart = (/1,kjstr/), kcount = (/jpiglo, kjcnt/) ) 
     1061         zbot(:,:) = zbot(:,:) * zbdy(:,:) 
     1062      ENDIF 
     1063      ! 
     1064      ldoce = zbot > 0. 
     1065      numout = inumsave 
     1066      ! 
     1067   END SUBROUTINE mpp_init_readbot_strip 
    6761068 
    6771069 
     
    7201112      ! 
    7211113   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 
    8341114 
    8351115 
     
    8961176   END SUBROUTINE mpp_init_nfdcom 
    8971177 
    898     
     1178 
    8991179#endif 
    9001180 
Note: See TracChangeset for help on using the changeset viewer.