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.

mppini.F90 in NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/LBC – NEMO

Context Navigation

source: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/LBC/mppini.F90 @ 10419

Last change on this file since 10419 was 10409, checked in by smasson, 6 years ago
dev_r10164_HPC09_ESIWACE_PREP_MERGE: action 6 add print when using non-optimal domain decomposition, see #2133
Property svn:keywords set to `Id`
File size: 53.9 KB

Rev	Line
[3]	1	MODULE mppini
[9019]	2	!!======================================================================
[3]	3	!! * MODULE mppini *
	4	!! Ocean initialization : distributed memory computing initialization
[9019]	5	!!======================================================================
	6	!! History : 6.0 ! 1994-11 (M. Guyon) Original code
	7	!! OPA 7.0 ! 1995-04 (J. Escobar, M. Imbard)
	8	!! 8.0 ! 1998-05 (M. Imbard, J. Escobar, L. Colombet ) SHMEM and MPI versions
	9	!! NEMO 1.0 ! 2004-01 (G. Madec, J.M Molines) F90 : free form , north fold jpni > 1
[9436]	10	!! 3.4 ! 2011-10 (A. C. Coward, NOCS & J. Donners, PRACE) add mpp_init_nfdcom
	11	!! 3. ! 2013-06 (I. Epicoco, S. Mocavero, CMCC) mpp_init_nfdcom: setup avoiding MPI communication
[9019]	12	!! 4.0 ! 2016-06 (G. Madec) use domain configuration file instead of bathymetry file
	13	!! 4.0 ! 2017-06 (J.M. Molines, T. Lovato) merge of mppini and mppini_2
	14	!!----------------------------------------------------------------------
[3]	15
	16	!!----------------------------------------------------------------------
[9436]	17	!! mpp_init : Lay out the global domain over processors with/without land processor elimination
	18	!! mpp_init_mask : Read global bathymetric information to facilitate land suppression
	19	!! mpp_init_ioipsl : IOIPSL initialization in mpp
	20	!! mpp_init_partition: Calculate MPP domain decomposition
	21	!! factorise : Calculate the factors of the no. of MPI processes
	22	!! mpp_init_nfdcom : Setup for north fold exchanges with explicit point-to-point messaging
[3]	23	!!----------------------------------------------------------------------
[9019]	24	USE dom_oce ! ocean space and time domain
	25	USE bdy_oce ! open BounDarY
	26	!
[9436]	27	USE lbcnfd , ONLY : isendto, nsndto, nfsloop, nfeloop ! Setup of north fold exchanges
[9019]	28	USE lib_mpp ! distribued memory computing library
	29	USE iom ! nemo I/O library
	30	USE ioipsl ! I/O IPSL library
	31	USE in_out_manager ! I/O Manager
[3]	32
	33	IMPLICIT NONE
	34	PRIVATE
	35
	36	PUBLIC mpp_init ! called by opa.F90
	37
[10330]	38	INTEGER :: numbot = -1 ! 'bottom_level' local logical unit
	39	INTEGER :: numbdy = -1 ! 'bdy_msk' local logical unit
	40
[3]	41	!!----------------------------------------------------------------------
[9598]	42	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[1152]	43	!! $Id$
[10068]	44	!! Software governed by the CeCILL license (see ./LICENSE)
[3]	45	!!----------------------------------------------------------------------
	46	CONTAINS
	47
[1556]	48	#if ! defined key_mpp_mpi
[3]	49	!!----------------------------------------------------------------------
	50	!! Default option : shared memory computing
	51	!!----------------------------------------------------------------------
	52
	53	SUBROUTINE mpp_init
	54	!!----------------------------------------------------------------------
	55	!! * ROUTINE mpp_init *
	56	!!
	57	!! ** Purpose : Lay out the global domain over processors.
	58	!!
	59	!! ** Method : Shared memory computing, set the local processor
[9019]	60	!! variables to the value of the global domain
[3]	61	!!----------------------------------------------------------------------
[9019]	62	!
[9436]	63	jpimax = jpiglo
	64	jpjmax = jpjglo
	65	jpi = jpiglo
	66	jpj = jpjglo
[9799]	67	jpk = jpkglo
[9436]	68	jpim1 = jpi-1 ! inner domain indices
	69	jpjm1 = jpj-1 ! " "
	70	jpkm1 = MAX( 1, jpk-1 ) ! " "
	71	jpij = jpi*jpj
	72	jpni = 1
	73	jpnj = 1
	74	jpnij = jpni*jpnj
[9019]	75	nimpp = 1 !
[3]	76	njmpp = 1
	77	nlci = jpi
	78	nlcj = jpj
	79	nldi = 1
	80	nldj = 1
	81	nlei = jpi
	82	nlej = jpj
	83	nbondi = 2
	84	nbondj = 2
[352]	85	nidom = FLIO_DOM_NONE
[1344]	86	npolj = jperio
[9667]	87	l_Iperio = jpni == 1 .AND. (jperio == 1 .OR. jperio == 4 .OR. jperio == 6 .OR. jperio == 7)
	88	l_Jperio = jpnj == 1 .AND. (jperio == 2 .OR. jperio == 7)
[9019]	89	!
[3]	90	IF(lwp) THEN
	91	WRITE(numout,*)
[9019]	92	WRITE(numout,*) 'mpp_init : NO massively parallel processing'
	93	WRITE(numout,*) '~~~~~~~~ '
[9667]	94	WRITE(numout,*) ' l_Iperio = ', l_Iperio, ' l_Jperio = ', l_Jperio
	95	WRITE(numout,*) ' npolj = ', npolj , ' njmpp = ', njmpp
[3]	96	ENDIF
[9019]	97	!
	98	IF( jpni /= 1 .OR. jpnj /= 1 .OR. jpnij /= 1 ) &
	99	CALL ctl_stop( 'mpp_init: equality jpni = jpnj = jpnij = 1 is not satisfied', &
	100	& 'the domain is lay out for distributed memory computing!' )
	101	!
[3]	102	END SUBROUTINE mpp_init
	103
	104	#else
	105	!!----------------------------------------------------------------------
[9019]	106	!! 'key_mpp_mpi' MPI massively parallel processing
[3]	107	!!----------------------------------------------------------------------
	108
[9449]	109
[3]	110	SUBROUTINE mpp_init
	111	!!----------------------------------------------------------------------
	112	!! * ROUTINE mpp_init *
	113	!!
	114	!! ** Purpose : Lay out the global domain over processors.
[9019]	115	!! If land processors are to be eliminated, this program requires the
	116	!! presence of the domain configuration file. Land processors elimination
	117	!! is performed if jpni x jpnj /= jpnij. In this case, using the MPP_PREP
	118	!! preprocessing tool, help for defining the best cutting out.
[3]	119	!!
	120	!! ** Method : Global domain is distributed in smaller local domains.
	121	!! Periodic condition is a function of the local domain position
	122	!! (global boundary or neighbouring domain) and of the global
	123	!! periodic
	124	!! Type : jperio global periodic condition
	125	!!
[9019]	126	!! ** Action : - set domain parameters
[3]	127	!! nimpp : longitudinal index
	128	!! njmpp : latitudinal index
	129	!! narea : number for local area
	130	!! nlci : first dimension
	131	!! nlcj : second dimension
	132	!! nbondi : mark for "east-west local boundary"
	133	!! nbondj : mark for "north-south local boundary"
	134	!! nproc : number for local processor
	135	!! noea : number for local neighboring processor
	136	!! nowe : number for local neighboring processor
	137	!! noso : number for local neighboring processor
	138	!! nono : number for local neighboring processor
	139	!!----------------------------------------------------------------------
[9019]	140	INTEGER :: ji, jj, jn, jproc, jarea ! dummy loop indices
[10345]	141	INTEGER :: inijmin
[10330]	142	INTEGER :: i2add
[9019]	143	INTEGER :: inum ! local logical unit
[10330]	144	INTEGER :: idir, ifreq, icont ! local integers
[9019]	145	INTEGER :: ii, il1, ili, imil ! - -
	146	INTEGER :: ij, il2, ilj, ijm1 ! - -
	147	INTEGER :: iino, ijno, iiso, ijso ! - -
	148	INTEGER :: iiea, ijea, iiwe, ijwe ! - -
[10330]	149	INTEGER :: iarea0 ! - -
	150	INTEGER :: ierr, ios !
[10409]	151	INTEGER :: inbi, inbj, iimax, ijmax !
[9444]	152	INTEGER, ALLOCATABLE, DIMENSION(:) :: iin, ii_nono, ii_noea ! 1D workspace
	153	INTEGER, ALLOCATABLE, DIMENSION(:) :: ijn, ii_noso, ii_nowe ! - -
	154	INTEGER, ALLOCATABLE, DIMENSION(:,:) :: iimppt, ilci, ibondi, ipproc ! 2D workspace
	155	INTEGER, ALLOCATABLE, DIMENSION(:,:) :: ijmppt, ilcj, ibondj, ipolj ! - -
	156	INTEGER, ALLOCATABLE, DIMENSION(:,:) :: ilei, ildi, iono, ioea ! - -
	157	INTEGER, ALLOCATABLE, DIMENSION(:,:) :: ilej, ildj, ioso, iowe ! - -
[10330]	158	LOGICAL, ALLOCATABLE, DIMENSION(:,:) :: llisoce ! - -
	159	NAMELIST/nambdy/ ln_bdy, nb_bdy, ln_coords_file, cn_coords_file, &
	160	& ln_mask_file, cn_mask_file, cn_dyn2d, nn_dyn2d_dta, &
	161	& cn_dyn3d, nn_dyn3d_dta, cn_tra, nn_tra_dta, &
	162	& ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp, rn_time_dmp_out, &
	163	& cn_ice, nn_ice_dta, &
	164	& rn_ice_tem, rn_ice_sal, rn_ice_age, &
	165	& ln_vol, nn_volctl, nn_rimwidth, nb_jpk_bdy
[3]	166	!!----------------------------------------------------------------------
[9436]	167
[10330]	168	! do we need to take into account bdy_msk?
	169	REWIND( numnam_ref ) ! Namelist nambdy in reference namelist : BDY
	170	READ ( numnam_ref, nambdy, IOSTAT = ios, ERR = 903)
	171	903 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy in reference namelist (mppini)', lwp )
	172	REWIND( numnam_cfg ) ! Namelist nambdy in configuration namelist : BDY
	173	READ ( numnam_cfg, nambdy, IOSTAT = ios, ERR = 904 )
	174	904 IF( ios > 0 ) CALL ctl_nam ( ios , 'nambdy in configuration namelist (mppini)', lwp )
	175	!
	176	IF( ln_read_cfg ) CALL iom_open( cn_domcfg, numbot )
	177	IF( ln_bdy .AND. ln_mask_file ) CALL iom_open( cn_mask_file, numbdy )
	178	!
	179	! 1. Dimension arrays for subdomains
	180	! -----------------------------------
	181	!
[9436]	182	! If dimensions of processor grid weren't specified in the namelist file
	183	! then we calculate them here now that we have our communicator size
[10409]	184	IF( jpni < 1 .OR. jpnj < 1 ) THEN
	185	CALL mpp_init_bestpartition( mppsize, jpni, jpnj )
	186	ELSE
	187	CALL mpp_init_bestpartition( mppsize, inbi, inbj )
	188	CALL mpp_basic_decomposition( jpni, jpnj, jpimax, jpjmax )
	189	CALL mpp_basic_decomposition( inbi, inbj, iimax, ijmax )
	190	IF( iimaxijmax < jpimaxjpjmax ) THEN
	191	WRITE(ctmp1,*) ' The chosen domain decomposition ', jpni, jpnj
	192	WRITE(ctmp2,) ' has larger MPI subdomains ', jpimax, jpjmax, jpimaxjpjmax
	193	WRITE(ctmp3,*) ' than the following domain decompostion ', inbi, inbj
	194	WRITE(ctmp4,) ' which MPI subdomains size is ', iimax, ijmax, iimaxijmax
	195	CALL ctl_warn( 'mpp_init:', '~~~~~~~~ ', ctmp1, ctmp2, ctmp3, ctmp4, ' --- YOU ARE WASTING CPU... ---' )
	196	ENDIF
	197	ENDIF
[10330]	198
	199	! look for land mpi subdomains...
	200	ALLOCATE( llisoce(jpni,jpnj) )
	201	CALL mpp_init_isoce( jpni, jpnj, llisoce )
	202	inijmin = COUNT( llisoce ) ! number of oce subdomains
	203
	204	IF( mppsize < inijmin ) THEN
	205	WRITE(ctmp1,*) ' With this specified domain decomposition: jpni =', jpni, ' jpnj =', jpnj
	206	WRITE(ctmp2,) ' we can eliminate only ', jpnijpnj - inijmin, ' land mpi subdomains therefore '
	207	WRITE(ctmp3,*) ' the number of ocean mpi subdomains (', inijmin,') exceed the number of MPI processes:', mppsize
	208	WRITE(ctmp4,*) ' ==>>> There is the list of best domain decompositions you should use: '
	209	CALL ctl_stop( 'mpp_init:', '~~~~~~~~ ', ctmp1, ctmp2, ctmp3, ctmp4 )
	210	CALL mpp_init_bestpartition( mppsize, ldlist = .TRUE. ) ! must be done by all core
	211	CALL ctl_stop( 'STOP' )
	212	ENDIF
	213
	214	IF( mppsize > jpni*jpnj ) THEN
	215	WRITE(ctmp1,*) ' The number of mpi processes: ', mppsize
	216	WRITE(ctmp2,) ' exceeds the maximum number of subdomains (ocean+land) = ', jpnijpnj
	217	WRITE(ctmp3,*) ' defined by the following domain decomposition: jpni =', jpni, ' jpnj =', jpnj
	218	WRITE(ctmp4,*) ' ==>>> There is the list of best domain decompositions you should use: '
	219	CALL ctl_stop( 'mpp_init:', '~~~~~~~~ ', ctmp1, ctmp2, ctmp3, ctmp4 )
	220	CALL mpp_init_bestpartition( mppsize, ldlist = .TRUE. ) ! must be done by all core
	221	CALL ctl_stop( 'STOP' )
	222	ENDIF
	223
	224	jpnij = mppsize ! force jpnij definition <-- remove as much land subdomains as needed to reach this condition
	225	IF( mppsize > inijmin ) THEN
	226	WRITE(ctmp1,*) ' The number of mpi processes: ', mppsize
	227	WRITE(ctmp2,*) ' exceeds the maximum number of ocean subdomains = ', inijmin
	228	WRITE(ctmp3,) ' we suppressed ', jpnijpnj - mppsize, ' land subdomains '
	229	WRITE(ctmp4,*) ' BUT we had to keep ', mppsize - inijmin, ' land subdomains that are useless...'
	230	CALL ctl_warn( 'mpp_init:', '~~~~~~~~ ', ctmp1, ctmp2, ctmp3, ctmp4 )
	231	ELSE ! mppsize = inijmin
	232	IF(lwp) THEN
	233	WRITE(numout,*) 'mpp_init: You use an optimal domain decomposition'
	234	WRITE(numout,*) '~~~~~~~~ '
	235	WRITE(numout,*) ' Number of mpi processes: ', mppsize
	236	WRITE(numout,*) ' Number of ocean subdomains = ', inijmin
	237	WRITE(numout,) ' Number of suppressed land subdomains = ', jpnijpnj - inijmin
	238	ENDIF
	239	ENDIF
	240
	241	IF( numbot /= -1 ) CALL iom_close( numbot )
	242	IF( numbdy /= -1 ) CALL iom_close( numbdy )
	243
[9436]	244	ALLOCATE( nfiimpp(jpni,jpnj), nfipproc(jpni,jpnj), nfilcit(jpni,jpnj) , &
	245	& nimppt(jpnij) , ibonit(jpnij) , nlcit(jpnij) , nlcjt(jpnij) , &
	246	& njmppt(jpnij) , ibonjt(jpnij) , nldit(jpnij) , nldjt(jpnij) , &
[9444]	247	& nleit(jpnij) , nlejt(jpnij) , &
	248	& iin(jpnij), ii_nono(jpnij), ii_noea(jpnij), &
	249	& ijn(jpnij), ii_noso(jpnij), ii_nowe(jpnij), &
	250	& iimppt(jpni,jpnj), ilci(jpni,jpnj), ibondi(jpni,jpnj), ipproc(jpni,jpnj), &
	251	& ijmppt(jpni,jpnj), ilcj(jpni,jpnj), ibondj(jpni,jpnj), ipolj(jpni,jpnj), &
	252	& ilei(jpni,jpnj), ildi(jpni,jpnj), iono(jpni,jpnj), ioea(jpni,jpnj), &
	253	& ilej(jpni,jpnj), ildj(jpni,jpnj), ioso(jpni,jpnj), iowe(jpni,jpnj), &
	254	& STAT=ierr )
[10297]	255	CALL mpp_sum( 'mppini', ierr )
[9436]	256	IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'mpp_init: unable to allocate standard ocean arrays' )
	257
[9508]	258	#if defined key_agrif
[9436]	259	IF( .NOT. Agrif_Root() ) THEN ! AGRIF children: specific setting (cf. agrif_user.F90)
[9449]	260	IF( jpiglo /= nbcellsx + 2 + 2*nbghostcells ) &
	261	CALL ctl_stop( 'STOP', 'mpp_init: Agrif children requires jpiglo == nbcellsx + 2 + 2*nbghostcells' )
	262	IF( jpjglo /= nbcellsy + 2 + 2*nbghostcells ) &
	263	CALL ctl_stop( 'STOP', 'mpp_init: Agrif children requires jpjglo == nbcellsy + 2 + 2*nbghostcells' )
	264	IF( ln_use_jattr ) CALL ctl_stop( 'STOP', 'mpp_init:Agrif children requires ln_use_jattr = .false. ' )
[9436]	265	ENDIF
[9508]	266	#endif
[9436]	267	!
[10330]	268	! 2. Index arrays for subdomains
[3]	269	! -----------------------------------
[9019]	270	!
	271	nreci = 2 * nn_hls
	272	nrecj = 2 * nn_hls
[10330]	273	CALL mpp_basic_decomposition( jpni, jpnj, jpimax, jpjmax, iimppt, ijmppt, ilci, ilcj )
	274	nfiimpp(:,:) = iimppt(:,:)
	275	nfilcit(:,:) = ilci(:,:)
[9019]	276	!
	277	IF(lwp) THEN
	278	WRITE(numout,*)
[10330]	279	WRITE(numout,*) 'MPI Message Passing MPI - domain lay out over processors'
	280	WRITE(numout,*)
[9019]	281	WRITE(numout,*) ' defines mpp subdomains'
[10330]	282	WRITE(numout,*) ' jpni = ', jpni
	283	WRITE(numout,*) ' jpnj = ', jpnj
[9019]	284	WRITE(numout,*)
[10330]	285	WRITE(numout,*) ' sum ilci(i,1) = ', sum(ilci(:,1)), ' jpiglo = ', jpiglo
	286	WRITE(numout,*) ' sum ilcj(1,j) = ', sum(ilcj(1,:)), ' jpjglo = ', jpjglo
[9019]	287	ENDIF
[10330]	288
[9019]	289	! 3. Subdomain description in the Regular Case
	290	! --------------------------------------------
[9667]	291	! specific cases where there is no communication -> must do the periodicity by itself
	292	! Warning: because of potential land-area suppression, do not use nbond[ij] == 2
	293	l_Iperio = jpni == 1 .AND. (jperio == 1 .OR. jperio == 4 .OR. jperio == 6 .OR. jperio == 7)
	294	l_Jperio = jpnj == 1 .AND. (jperio == 2 .OR. jperio == 7)
	295
[9019]	296	DO jarea = 1, jpni*jpnj
[10330]	297	!
[9564]	298	iarea0 = jarea - 1
	299	ii = 1 + MOD(iarea0,jpni)
	300	ij = 1 + iarea0/jpni
[9019]	301	ili = ilci(ii,ij)
	302	ilj = ilcj(ii,ij)
[9667]	303	ibondi(ii,ij) = 0 ! default: has e-w neighbours
	304	IF( ii == 1 ) ibondi(ii,ij) = -1 ! first column, has only e neighbour
	305	IF( ii == jpni ) ibondi(ii,ij) = 1 ! last column, has only w neighbour
	306	IF( jpni == 1 ) ibondi(ii,ij) = 2 ! has no e-w neighbour
	307	ibondj(ii,ij) = 0 ! default: has n-s neighbours
	308	IF( ij == 1 ) ibondj(ii,ij) = -1 ! first row, has only n neighbour
	309	IF( ij == jpnj ) ibondj(ii,ij) = 1 ! last row, has only s neighbour
	310	IF( jpnj == 1 ) ibondj(ii,ij) = 2 ! has no n-s neighbour
[9019]	311
[9436]	312	! Subdomain neighbors (get their zone number): default definition
[9564]	313	ioso(ii,ij) = iarea0 - jpni
	314	iowe(ii,ij) = iarea0 - 1
	315	ioea(ii,ij) = iarea0 + 1
	316	iono(ii,ij) = iarea0 + jpni
[9019]	317	ildi(ii,ij) = 1 + nn_hls
	318	ilei(ii,ij) = ili - nn_hls
	319	ildj(ii,ij) = 1 + nn_hls
	320	ilej(ii,ij) = ilj - nn_hls
	321
[9436]	322	! East-West periodicity: change ibondi, ioea, iowe
[9069]	323	IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 .OR. jperio == 7 ) THEN
[9667]	324	IF( jpni /= 1 ) ibondi(ii,ij) = 0 ! redefine: all have e-w neighbours
	325	IF( ii == 1 ) iowe(ii,ij) = iarea0 + (jpni-1) ! redefine: first column, address of w neighbour
	326	IF( ii == jpni ) ioea(ii,ij) = iarea0 - (jpni-1) ! redefine: last column, address of e neighbour
[9019]	327	ENDIF
[9436]	328
[9667]	329	! Simple North-South periodicity: change ibondj, ioso, iono
	330	IF( jperio == 2 .OR. jperio == 7 ) THEN
	331	IF( jpnj /= 1 ) ibondj(ii,ij) = 0 ! redefine: all have n-s neighbours
	332	IF( ij == 1 ) ioso(ii,ij) = iarea0 + jpni * (jpnj-1) ! redefine: first row, address of s neighbour
	333	IF( ij == jpnj ) iono(ii,ij) = iarea0 - jpni * (jpnj-1) ! redefine: last row, address of n neighbour
	334	ENDIF
	335
[9436]	336	! North fold: define ipolj, change iono. Warning: we do not change ibondj...
[9019]	337	ipolj(ii,ij) = 0
	338	IF( jperio == 3 .OR. jperio == 4 ) THEN
	339	ijm1 = jpni*(jpnj-1)
	340	imil = ijm1+(jpni+1)/2
	341	IF( jarea > ijm1 ) ipolj(ii,ij) = 3
	342	IF( MOD(jpni,2) == 1 .AND. jarea == imil ) ipolj(ii,ij) = 4
	343	IF( ipolj(ii,ij) == 3 ) iono(ii,ij) = jpni*jpnj-jarea+ijm1 ! MPI rank of northern neighbour
	344	ENDIF
	345	IF( jperio == 5 .OR. jperio == 6 ) THEN
	346	ijm1 = jpni*(jpnj-1)
	347	imil = ijm1+(jpni+1)/2
	348	IF( jarea > ijm1) ipolj(ii,ij) = 5
	349	IF( MOD(jpni,2) == 1 .AND. jarea == imil ) ipolj(ii,ij) = 6
	350	IF( ipolj(ii,ij) == 5) iono(ii,ij) = jpni*jpnj-jarea+ijm1 ! MPI rank of northern neighbour
	351	ENDIF
	352	!
[10330]	353	END DO
	354
	355	! 4. deal with land subdomains
	356	! ----------------------------
	357	!
	358	! specify which subdomains are oce subdomains; other are land subdomains
	359	ipproc(:,:) = -1
	360	icont = -1
	361	DO jarea = 1, jpni*jpnj
	362	iarea0 = jarea - 1
	363	ii = 1 + MOD(iarea0,jpni)
	364	ij = 1 + iarea0/jpni
	365	IF( llisoce(ii,ij) ) THEN
[9019]	366	icont = icont + 1
	367	ipproc(ii,ij) = icont
	368	iin(icont+1) = ii
	369	ijn(icont+1) = ij
	370	ENDIF
[3]	371	END DO
[10330]	372	! if needed add some land subdomains to reach jpnij active subdomains
	373	i2add = jpnij - inijmin
	374	DO jarea = 1, jpni*jpnj
	375	iarea0 = jarea - 1
	376	ii = 1 + MOD(iarea0,jpni)
	377	ij = 1 + iarea0/jpni
	378	IF( .NOT. llisoce(ii,ij) .AND. i2add > 0 ) THEN
	379	icont = icont + 1
	380	ipproc(ii,ij) = icont
	381	iin(icont+1) = ii
	382	ijn(icont+1) = ij
	383	i2add = i2add - 1
	384	ENDIF
	385	END DO
[9019]	386	nfipproc(:,:) = ipproc(:,:)
[6412]	387
[10330]	388	! neighbour treatment: change ibondi, ibondj if next to a land zone
[9019]	389	DO jarea = 1, jpni*jpnj
	390	ii = 1 + MOD( jarea-1 , jpni )
	391	ij = 1 + (jarea-1) / jpni
[9667]	392	! land-only area with an active n neigbour
[9019]	393	IF ( ipproc(ii,ij) == -1 .AND. 0 <= iono(ii,ij) .AND. iono(ii,ij) <= jpni*jpnj-1 ) THEN
[9667]	394	iino = 1 + MOD( iono(ii,ij) , jpni ) ! ii index of this n neigbour
	395	ijno = 1 + iono(ii,ij) / jpni ! ij index of this n neigbour
	396	! In case of north fold exchange: I am the n neigbour of my n neigbour!! (#1057)
	397	! --> for northern neighbours of northern row processors (in case of north-fold)
	398	! need to reverse the LOGICAL direction of communication
	399	idir = 1 ! we are indeed the s neigbour of this n neigbour
	400	IF( ij == jpnj .AND. ijno == jpnj ) idir = -1 ! both are on the last row, we are in fact the n neigbour
	401	IF( ibondj(iino,ijno) == idir ) ibondj(iino,ijno) = 2 ! this n neigbour had only a s/n neigbour -> no more
	402	IF( ibondj(iino,ijno) == 0 ) ibondj(iino,ijno) = -idir ! this n neigbour had both, n-s neighbours -> keep 1
[9019]	403	ENDIF
[9667]	404	! land-only area with an active s neigbour
[9019]	405	IF( ipproc(ii,ij) == -1 .AND. 0 <= ioso(ii,ij) .AND. ioso(ii,ij) <= jpni*jpnj-1 ) THEN
[9667]	406	iiso = 1 + MOD( ioso(ii,ij) , jpni ) ! ii index of this s neigbour
	407	ijso = 1 + ioso(ii,ij) / jpni ! ij index of this s neigbour
	408	IF( ibondj(iiso,ijso) == -1 ) ibondj(iiso,ijso) = 2 ! this s neigbour had only a n neigbour -> no more neigbour
	409	IF( ibondj(iiso,ijso) == 0 ) ibondj(iiso,ijso) = 1 ! this s neigbour had both, n-s neighbours -> keep s neigbour
[9019]	410	ENDIF
[9667]	411	! land-only area with an active e neigbour
[9019]	412	IF( ipproc(ii,ij) == -1 .AND. 0 <= ioea(ii,ij) .AND. ioea(ii,ij) <= jpni*jpnj-1 ) THEN
[9667]	413	iiea = 1 + MOD( ioea(ii,ij) , jpni ) ! ii index of this e neigbour
	414	ijea = 1 + ioea(ii,ij) / jpni ! ij index of this e neigbour
	415	IF( ibondi(iiea,ijea) == 1 ) ibondi(iiea,ijea) = 2 ! this e neigbour had only a w neigbour -> no more neigbour
	416	IF( ibondi(iiea,ijea) == 0 ) ibondi(iiea,ijea) = -1 ! this e neigbour had both, e-w neighbours -> keep e neigbour
[9019]	417	ENDIF
[9667]	418	! land-only area with an active w neigbour
[9019]	419	IF( ipproc(ii,ij) == -1 .AND. 0 <= iowe(ii,ij) .AND. iowe(ii,ij) <= jpni*jpnj-1) THEN
[9667]	420	iiwe = 1 + MOD( iowe(ii,ij) , jpni ) ! ii index of this w neigbour
	421	ijwe = 1 + iowe(ii,ij) / jpni ! ij index of this w neigbour
	422	IF( ibondi(iiwe,ijwe) == -1 ) ibondi(iiwe,ijwe) = 2 ! this w neigbour had only a e neigbour -> no more neigbour
	423	IF( ibondi(iiwe,ijwe) == 0 ) ibondi(iiwe,ijwe) = 1 ! this w neigbour had both, e-w neighbours -> keep w neigbour
[9019]	424	ENDIF
	425	END DO
[3]	426
[9436]	427	! Update il[de][ij] according to modified ibond[ij]
	428	! ----------------------
[9564]	429	DO jproc = 1, jpnij
	430	ii = iin(jproc)
	431	ij = ijn(jproc)
[9436]	432	IF( ibondi(ii,ij) == -1 .OR. ibondi(ii,ij) == 2 ) ildi(ii,ij) = 1
	433	IF( ibondi(ii,ij) == 1 .OR. ibondi(ii,ij) == 2 ) ilei(ii,ij) = ilci(ii,ij)
	434	IF( ibondj(ii,ij) == -1 .OR. ibondj(ii,ij) == 2 ) ildj(ii,ij) = 1
	435	IF( ibondj(ii,ij) == 1 .OR. ibondj(ii,ij) == 2 ) ilej(ii,ij) = ilcj(ii,ij)
	436	END DO
[10330]	437
	438	! 5. Subdomain print
	439	! ------------------
	440	IF(lwp) THEN
	441	ifreq = 4
	442	il1 = 1
	443	DO jn = 1, (jpni-1)/ifreq+1
	444	il2 = MIN(jpni,il1+ifreq-1)
	445	WRITE(numout,*)
	446	WRITE(numout,9400) ('***',ji=il1,il2-1)
	447	DO jj = jpnj, 1, -1
	448	WRITE(numout,9403) (' ',ji=il1,il2-1)
	449	WRITE(numout,9402) jj, (ilci(ji,jj),ilcj(ji,jj),ji=il1,il2)
	450	WRITE(numout,9404) (ipproc(ji,jj),ji=il1,il2)
	451	WRITE(numout,9403) (' ',ji=il1,il2-1)
	452	WRITE(numout,9400) ('***',ji=il1,il2-1)
	453	END DO
	454	WRITE(numout,9401) (ji,ji=il1,il2)
	455	il1 = il1+ifreq
	456	END DO
	457	9400 FORMAT(' *' ,20('***********',a3) )
	458	9403 FORMAT(' * ',20(' * ',a3) )
	459	9401 FORMAT(' ' ,20(' ',i3,' ') )
	460	9402 FORMAT(' ',i3,' * ',20(i3,' x',i3,' * ') )
	461	9404 FORMAT(' * ' ,20(' ',i3,' * ') )
	462	ENDIF
[9436]	463
[9019]	464	! just to save nono etc for all proc
[9436]	465	! warning ii*ij (zone) /= nproc (processors)!
	466	! ioso = zone number, ii_noso = proc number
[9019]	467	ii_noso(:) = -1
	468	ii_nono(:) = -1
	469	ii_noea(:) = -1
	470	ii_nowe(:) = -1
[9449]	471	DO jproc = 1, jpnij
	472	ii = iin(jproc)
	473	ij = ijn(jproc)
[9019]	474	IF( 0 <= ioso(ii,ij) .AND. ioso(ii,ij) <= (jpni*jpnj-1) ) THEN
	475	iiso = 1 + MOD( ioso(ii,ij) , jpni )
	476	ijso = 1 + ioso(ii,ij) / jpni
[9449]	477	ii_noso(jproc) = ipproc(iiso,ijso)
[9019]	478	ENDIF
	479	IF( 0 <= iowe(ii,ij) .AND. iowe(ii,ij) <= (jpni*jpnj-1) ) THEN
	480	iiwe = 1 + MOD( iowe(ii,ij) , jpni )
	481	ijwe = 1 + iowe(ii,ij) / jpni
[9449]	482	ii_nowe(jproc) = ipproc(iiwe,ijwe)
[9019]	483	ENDIF
	484	IF( 0 <= ioea(ii,ij) .AND. ioea(ii,ij) <= (jpni*jpnj-1) ) THEN
	485	iiea = 1 + MOD( ioea(ii,ij) , jpni )
	486	ijea = 1 + ioea(ii,ij) / jpni
[9449]	487	ii_noea(jproc)= ipproc(iiea,ijea)
[9019]	488	ENDIF
	489	IF( 0 <= iono(ii,ij) .AND. iono(ii,ij) <= (jpni*jpnj-1) ) THEN
	490	iino = 1 + MOD( iono(ii,ij) , jpni )
	491	ijno = 1 + iono(ii,ij) / jpni
[9449]	492	ii_nono(jproc)= ipproc(iino,ijno)
[9019]	493	ENDIF
	494	END DO
	495
	496	! 6. Change processor name
	497	! ------------------------
	498	ii = iin(narea)
	499	ij = ijn(narea)
	500	!
	501	! set default neighbours
	502	noso = ii_noso(narea)
	503	nowe = ii_nowe(narea)
	504	noea = ii_noea(narea)
	505	nono = ii_nono(narea)
	506	nlci = ilci(ii,ij)
	507	nldi = ildi(ii,ij)
	508	nlei = ilei(ii,ij)
[9436]	509	nlcj = ilcj(ii,ij)
[9019]	510	nldj = ildj(ii,ij)
	511	nlej = ilej(ii,ij)
	512	nbondi = ibondi(ii,ij)
	513	nbondj = ibondj(ii,ij)
	514	nimpp = iimppt(ii,ij)
[9436]	515	njmpp = ijmppt(ii,ij)
	516	jpi = nlci
	517	jpj = nlcj
	518	jpk = jpkglo ! third dim
	519	#if defined key_agrif
	520	! simple trick to use same vertical grid as parent but different number of levels:
	521	! Save maximum number of levels in jpkglo, then define all vertical grids with this number.
	522	! Suppress once vertical online interpolation is ok
[9449]	523	!!$ IF(.NOT.Agrif_Root()) jpkglo = Agrif_Parent( jpkglo )
[9436]	524	#endif
	525	jpim1 = jpi-1 ! inner domain indices
	526	jpjm1 = jpj-1 ! " "
	527	jpkm1 = MAX( 1, jpk-1 ) ! " "
	528	jpij = jpi*jpj ! jpi x j
[9019]	529	DO jproc = 1, jpnij
	530	ii = iin(jproc)
	531	ij = ijn(jproc)
	532	nlcit(jproc) = ilci(ii,ij)
	533	nldit(jproc) = ildi(ii,ij)
	534	nleit(jproc) = ilei(ii,ij)
[9436]	535	nlcjt(jproc) = ilcj(ii,ij)
[9019]	536	nldjt(jproc) = ildj(ii,ij)
	537	nlejt(jproc) = ilej(ii,ij)
[9436]	538	ibonit(jproc) = ibondi(ii,ij)
	539	ibonjt(jproc) = ibondj(ii,ij)
	540	nimppt(jproc) = iimppt(ii,ij)
	541	njmppt(jproc) = ijmppt(ii,ij)
[9019]	542	END DO
[3]	543
[9019]	544	! Save processor layout in ascii file
	545	IF (lwp) THEN
[7646]	546	CALL ctl_opn( inum, 'layout.dat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE., narea )
[9019]	547	WRITE(inum,'(a)') ' jpnij jpimax jpjmax jpk jpiglo jpjglo'//&
[9436]	548	& ' ( local: narea jpi jpj )'
[9019]	549	WRITE(inum,'(6i8,a,3i8,a)') jpnij,jpimax,jpjmax,jpk,jpiglo,jpjglo,&
	550	& ' ( local: ',narea,jpi,jpj,' )'
[9436]	551	WRITE(inum,'(a)') 'nproc nlci nlcj nldi nldj nlei nlej nimp njmp nono noso nowe noea nbondi nbondj '
[9019]	552
	553	DO jproc = 1, jpnij
[9436]	554	WRITE(inum,'(13i5,2i7)') jproc-1, nlcit (jproc), nlcjt (jproc), &
	555	& nldit (jproc), nldjt (jproc), &
	556	& nleit (jproc), nlejt (jproc), &
	557	& nimppt (jproc), njmppt (jproc), &
	558	& ii_nono(jproc), ii_noso(jproc), &
	559	& ii_nowe(jproc), ii_noea(jproc), &
[9446]	560	& ibonit (jproc), ibonjt (jproc)
[7646]	561	END DO
[3]	562	END IF
	563
[9019]	564	! ! north fold parameter
	565	! Defined npolj, either 0, 3 , 4 , 5 , 6
	566	! In this case the important thing is that npolj /= 0
	567	! Because if we go through these line it is because jpni >1 and thus
	568	! we must use lbcnorthmpp, which tests only npolj =0 or npolj /= 0
[3]	569	npolj = 0
[9019]	570	ij = ijn(narea)
[3]	571	IF( jperio == 3 .OR. jperio == 4 ) THEN
[9019]	572	IF( ij == jpnj ) npolj = 3
[3]	573	ENDIF
	574	IF( jperio == 5 .OR. jperio == 6 ) THEN
[9019]	575	IF( ij == jpnj ) npolj = 5
[3]	576	ENDIF
[9019]	577	!
[9436]	578	nproc = narea-1
[3]	579	IF(lwp) THEN
[6412]	580	WRITE(numout,*)
[9169]	581	WRITE(numout,*) ' resulting internal parameters : '
	582	WRITE(numout,*) ' nproc = ', nproc
	583	WRITE(numout,*) ' nowe = ', nowe , ' noea = ', noea
	584	WRITE(numout,*) ' nono = ', nono , ' noso = ', noso
	585	WRITE(numout,*) ' nbondi = ', nbondi
	586	WRITE(numout,*) ' nbondj = ', nbondj
	587	WRITE(numout,*) ' npolj = ', npolj
[9667]	588	WRITE(numout,*) ' l_Iperio = ', l_Iperio
	589	WRITE(numout,*) ' l_Jperio = ', l_Jperio
[9169]	590	WRITE(numout,*) ' nlci = ', nlci
	591	WRITE(numout,*) ' nlcj = ', nlcj
	592	WRITE(numout,*) ' nimpp = ', nimpp
	593	WRITE(numout,*) ' njmpp = ', njmpp
	594	WRITE(numout,*) ' nreci = ', nreci
	595	WRITE(numout,*) ' nrecj = ', nrecj
	596	WRITE(numout,*) ' nn_hls = ', nn_hls
[3]	597	ENDIF
	598
[9019]	599	! ! Prepare mpp north fold
[6412]	600	IF( jperio >= 3 .AND. jperio <= 6 .AND. jpni > 1 ) THEN
[3]	601	CALL mpp_ini_north
[9917]	602	IF (lwp) THEN
	603	WRITE(numout,*)
	604	WRITE(numout,*) ' ==>>> North fold boundary prepared for jpni >1'
	605	! additional prints in layout.dat
	606	WRITE(inum,*)
	607	WRITE(inum,*)
	608	WRITE(inum,*) 'number of subdomains located along the north fold : ', ndim_rank_north
	609	WRITE(inum,*) 'Rank of the subdomains located along the north fold : ', ndim_rank_north
	610	DO jproc = 1, ndim_rank_north, 5
	611	WRITE(inum,*) nrank_north( jproc:MINVAL( (/jproc+4,ndim_rank_north/) ) )
	612	END DO
	613	ENDIF
[6412]	614	ENDIF
[9019]	615	!
	616	CALL mpp_init_ioipsl ! Prepare NetCDF output file (if necessary)
	617	!
[9917]	618	IF( ln_nnogather ) THEN
	619	CALL mpp_init_nfdcom ! northfold neighbour lists
	620	IF (lwp) THEN
	621	WRITE(inum,*)
	622	WRITE(inum,*)
	623	WRITE(inum,*) 'north fold exchanges with explicit point-to-point messaging :'
	624	WRITE(inum,*) 'nfsloop : ', nfsloop
	625	WRITE(inum,*) 'nfeloop : ', nfeloop
	626	WRITE(inum,*) 'nsndto : ', nsndto
	627	WRITE(inum,*) 'isendto : ', isendto
	628	ENDIF
	629	ENDIF
[9436]	630	!
[9917]	631	IF (lwp) CLOSE(inum)
	632	!
[9444]	633	DEALLOCATE(iin, ijn, ii_nono, ii_noea, ii_noso, ii_nowe, &
	634	& iimppt, ijmppt, ibondi, ibondj, ipproc, ipolj, &
	635	& ilci, ilcj, ilei, ilej, ildi, ildj, &
[10330]	636	& iono, ioea, ioso, iowe, llisoce)
[9444]	637	!
[9019]	638	END SUBROUTINE mpp_init
[3]	639
	640
[10330]	641	SUBROUTINE mpp_basic_decomposition( knbi, knbj, kimax, kjmax, kimppt, kjmppt, klci, klcj)
[9019]	642	!!----------------------------------------------------------------------
[10330]	643	!! * ROUTINE mpp_basic_decomposition *
	644	!!
	645	!! ** Purpose : Lay out the global domain over processors.
[9019]	646	!!
[10330]	647	!! ** Method : Global domain is distributed in smaller local domains.
[9019]	648	!!
[10330]	649	!! ** Action : - set for all knbi*knbj domains:
	650	!! kimppt : longitudinal index
	651	!! kjmppt : latitudinal index
	652	!! klci : first dimension
	653	!! klcj : second dimension
[9019]	654	!!----------------------------------------------------------------------
[10409]	655	INTEGER, INTENT(in ) :: knbi, knbj
	656	INTEGER, INTENT( out) :: kimax, kjmax
	657	INTEGER, DIMENSION(knbi,knbj), OPTIONAL, INTENT( out) :: kimppt, kjmppt
	658	INTEGER, DIMENSION(knbi,knbj), OPTIONAL, INTENT( out) :: klci, klcj
[10330]	659	!
	660	INTEGER :: ji, jj
	661	INTEGER :: iresti, irestj
	662	INTEGER :: ireci, irecj
[9019]	663	!!----------------------------------------------------------------------
	664	!
[10330]	665	#if defined key_nemocice_decomp
[10409]	666	kimax = ( nx_global+2-2nn_hls + (knbi-1) ) / knbi + 2nn_hls ! first dim.
	667	kjmax = ( ny_global+2-2nn_hls + (knbj-1) ) / knbj + 2nn_hls ! second dim.
[10330]	668	#else
	669	kimax = ( jpiglo - 2nn_hls + (knbi-1) ) / knbi + 2nn_hls ! first dim.
	670	kjmax = ( jpjglo - 2nn_hls + (knbj-1) ) / knbj + 2nn_hls ! second dim.
	671	#endif
[10409]	672	IF( .NOT. PRESENT(kimppt) ) RETURN
[9019]	673	!
[10330]	674	! 1. Dimension arrays for subdomains
	675	! -----------------------------------
	676	! Computation of local domain sizes klci() klcj()
	677	! These dimensions depend on global sizes knbi,knbj and jpiglo,jpjglo
	678	! The subdomains are squares lesser than or equal to the global
	679	! dimensions divided by the number of processors minus the overlap array.
[9019]	680	!
[10330]	681	ireci = 2 * nn_hls
	682	irecj = 2 * nn_hls
	683	iresti = 1 + MOD( jpiglo - ireci -1 , knbi )
	684	irestj = 1 + MOD( jpjglo - irecj -1 , knbj )
[9168]	685	!
[10330]	686	! Need to use kimax and kjmax here since jpi and jpj not yet defined
	687	#if defined key_nemocice_decomp
	688	! Change padding to be consistent with CICE
	689	klci(1:knbi-1 ,:) = kimax
	690	klci(knbi ,:) = jpiglo - (knbi - 1) * (kimax - nreci)
	691	klcj(:, 1:knbj-1) = kjmax
	692	klcj(:, knbj) = jpjglo - (knbj - 1) * (kjmax - nrecj)
	693	#else
	694	klci(1:iresti ,:) = kimax
	695	klci(iresti+1:knbi ,:) = kimax-1
	696	klcj(:, 1:irestj) = kjmax
	697	klcj(:, irestj+1:knbj) = kjmax-1
	698	#endif
[3]	699
[10330]	700	! 2. Index arrays for subdomains
	701	! -------------------------------
	702	kimppt(:,:) = 1
	703	kjmppt(:,:) = 1
	704	!
	705	IF( knbi > 1 ) THEN
	706	DO jj = 1, knbj
	707	DO ji = 2, knbi
	708	kimppt(ji,jj) = kimppt(ji-1,jj) + klci(ji-1,jj) - ireci
	709	END DO
	710	END DO
[9019]	711	ENDIF
	712	!
[10330]	713	IF( knbj > 1 )THEN
	714	DO jj = 2, knbj
	715	DO ji = 1, knbi
	716	kjmppt(ji,jj) = kjmppt(ji,jj-1) + klcj(ji,jj-1) - irecj
	717	END DO
	718	END DO
	719	ENDIF
	720
	721	END SUBROUTINE mpp_basic_decomposition
[9019]	722
	723
[10330]	724	SUBROUTINE mpp_init_bestpartition( knbij, knbi, knbj, ldlist )
	725	!!----------------------------------------------------------------------
	726	!! * ROUTINE mpp_init_bestpartition *
	727	!!
	728	!! ** Purpose :
	729	!!
	730	!! ** Method :
	731	!!----------------------------------------------------------------------
	732	INTEGER, INTENT(in ) :: knbij ! total number if subdomains (knbi*knbj)
	733	INTEGER, OPTIONAL, INTENT( out) :: knbi, knbj ! number if subdomains along i and j (knbi and knbj)
	734	LOGICAL, OPTIONAL, INTENT(in ) :: ldlist ! .true.: print the list the best domain decompositions (with land)
	735	!
	736	INTEGER :: ji, jj, ii, iitarget
	737	INTEGER :: iszitst, iszjtst
	738	INTEGER :: isziref, iszjref
	739	INTEGER :: inbij, iszij
	740	INTEGER :: inbimax, inbjmax, inbijmax
	741	INTEGER :: isz0, isz1
	742	INTEGER, DIMENSION( :), ALLOCATABLE :: indexok
	743	INTEGER, DIMENSION( :), ALLOCATABLE :: inbi0, inbj0, inbij0 ! number of subdomains along i,j
	744	INTEGER, DIMENSION( :), ALLOCATABLE :: iszi0, iszj0, iszij0 ! max size of the subdomains along i,j
	745	INTEGER, DIMENSION( :), ALLOCATABLE :: inbi1, inbj1, inbij1 ! number of subdomains along i,j
	746	INTEGER, DIMENSION( :), ALLOCATABLE :: iszi1, iszj1, iszij1 ! max size of the subdomains along i,j
	747	LOGICAL :: llist
	748	LOGICAL, DIMENSION(:,:), ALLOCATABLE :: llmsk2d ! max size of the subdomains along i,j
	749	LOGICAL, DIMENSION(:,:), ALLOCATABLE :: llisoce ! - -
	750	REAL(wp):: zpropland
	751	!!----------------------------------------------------------------------
	752	!
	753	llist = .FALSE.
	754	IF( PRESENT(ldlist) ) llist = ldlist
	755
	756	CALL mpp_init_landprop( zpropland ) ! get the proportion of land point over the gloal domain
	757	inbij = NINT( REAL(knbij, wp) / ( 1.0 - zpropland ) ) ! define the largest possible value for jpni*jpnj
	758	!
	759	IF( llist ) THEN ; inbijmax = inbij*2
	760	ELSE ; inbijmax = inbij
	761	ENDIF
	762	!
	763	ALLOCATE(inbi0(inbijmax),inbj0(inbijmax),iszi0(inbijmax),iszj0(inbijmax))
	764	!
	765	inbimax = 0
	766	inbjmax = 0
	767	isziref = jpiglo*jpjglo+1
	768	iszjref = jpiglo*jpjglo+1
	769	!
	770	! get the list of knbi that gives a smaller jpimax than knbi-1
	771	! get the list of knbj that gives a smaller jpjmax than knbj-1
	772	DO ji = 1, inbijmax
	773	#if defined key_nemocice_decomp
	774	iszitst = ( nx_global+2-2nn_hls + (ji-1) ) / ji + 2nn_hls ! first dim.
	775	#else
	776	iszitst = ( jpiglo - 2nn_hls + (ji-1) ) / ji + 2nn_hls
	777	#endif
	778	IF( iszitst < isziref ) THEN
	779	isziref = iszitst
	780	inbimax = inbimax + 1
	781	inbi0(inbimax) = ji
	782	iszi0(inbimax) = isziref
	783	ENDIF
	784	#if defined key_nemocice_decomp
	785	iszjtst = ( ny_global+2-2nn_hls + (ji-1) ) / ji + 2nn_hls ! first dim.
	786	#else
	787	iszjtst = ( jpjglo - 2nn_hls + (ji-1) ) / ji + 2nn_hls
	788	#endif
	789	IF( iszjtst < iszjref ) THEN
	790	iszjref = iszjtst
	791	inbjmax = inbjmax + 1
	792	inbj0(inbjmax) = ji
	793	iszj0(inbjmax) = iszjref
	794	ENDIF
	795	END DO
	796
	797	! combine these 2 lists to get all possible knbi*knbj < inbijmax
	798	ALLOCATE( llmsk2d(inbimax,inbjmax) )
	799	DO jj = 1, inbjmax
	800	DO ji = 1, inbimax
	801	IF ( inbi0(ji) * inbj0(jj) <= inbijmax ) THEN ; llmsk2d(ji,jj) = .TRUE.
	802	ELSE ; llmsk2d(ji,jj) = .FALSE.
	803	ENDIF
	804	END DO
	805	END DO
	806	isz1 = COUNT(llmsk2d)
	807	ALLOCATE( inbi1(isz1), inbj1(isz1), iszi1(isz1), iszj1(isz1) )
	808	ii = 0
	809	DO jj = 1, inbjmax
	810	DO ji = 1, inbimax
[10357]	811	IF( llmsk2d(ji,jj) .EQV. .TRUE. ) THEN
[10330]	812	ii = ii + 1
	813	inbi1(ii) = inbi0(ji)
	814	inbj1(ii) = inbj0(jj)
	815	iszi1(ii) = iszi0(ji)
	816	iszj1(ii) = iszj0(jj)
	817	END IF
	818	END DO
	819	END DO
	820	DEALLOCATE( inbi0, inbj0, iszi0, iszj0 )
	821	DEALLOCATE( llmsk2d )
	822
	823	ALLOCATE( inbij1(isz1), iszij1(isz1) )
	824	inbij1(:) = inbi1(:) * inbj1(:)
	825	iszij1(:) = iszi1(:) * iszj1(:)
	826
	827	! if therr is no land and no print
	828	IF( .NOT. llist .AND. numbot == -1 .AND. numbdy == -1 ) THEN
	829	! get the smaller partition which gives the smallest subdomain size
	830	ii = MINLOC(inbij1, mask = iszij1 == MINVAL(iszij1), dim = 1)
	831	knbi = inbi1(ii)
	832	knbj = inbj1(ii)
	833	DEALLOCATE( inbi1, inbj1, inbij1, iszi1, iszj1, iszij1 )
	834	RETURN
	835	ENDIF
	836
	837	! extract only the partitions which reduce the subdomain size in comparison with smaller partitions
	838	ALLOCATE( indexok(isz1) ) ! to store indices of the best partitions
	839	isz0 = 0 ! number of best partitions
	840	inbij = 1 ! start with the min value of inbij1 => 1
	841	iszij = jpiglo*jpjglo+1 ! default: larger than global domain
	842	DO WHILE( inbij <= inbijmax ) ! if we did not reach the max of inbij1
	843	ii = MINLOC(iszij1, mask = inbij1 == inbij, dim = 1) ! warning: send back the first occurence if multiple results
	844	IF ( iszij1(ii) < iszij ) THEN
	845	isz0 = isz0 + 1
	846	indexok(isz0) = ii
	847	iszij = iszij1(ii)
	848	ENDIF
	849	inbij = MINVAL(inbij1, mask = inbij1 > inbij) ! warning: return largest integer value if mask = .false. everywhere
	850	END DO
	851	DEALLOCATE( inbij1, iszij1 )
	852
	853	! keep only the best partitions (sorted by increasing order of subdomains number and decreassing subdomain size)
	854	ALLOCATE( inbi0(isz0), inbj0(isz0), iszi0(isz0), iszj0(isz0) )
	855	DO ji = 1, isz0
	856	ii = indexok(ji)
	857	inbi0(ji) = inbi1(ii)
	858	inbj0(ji) = inbj1(ii)
	859	iszi0(ji) = iszi1(ii)
	860	iszj0(ji) = iszj1(ii)
	861	END DO
	862	DEALLOCATE( indexok, inbi1, inbj1, iszi1, iszj1 )
	863
	864	IF( llist ) THEN ! we print about 21 best partitions
	865	IF(lwp) THEN
	866	WRITE(numout,*)
	867	WRITE(numout, *) ' For your information:'
	868	WRITE(numout,'(a,i5,a)') ' list of the best partitions around ', knbij, ' mpi processes'
	869	WRITE(numout, *) ' --------------------------------------', '-----', '--------------'
	870	WRITE(numout,*)
	871	END IF
	872	iitarget = MINLOC( inbi0(:)inbj0(:), mask = inbi0(:)inbj0(:) >= knbij, dim = 1 )
	873	DO ji = MAX(1,iitarget-10), MIN(isz0,iitarget+10)
	874	ALLOCATE( llisoce(inbi0(ji), inbj0(ji)) )
	875	CALL mpp_init_isoce( inbi0(ji), inbj0(ji), llisoce ) ! Warning: must be call by all cores (call mpp_sum)
	876	inbij = COUNT(llisoce)
	877	DEALLOCATE( llisoce )
	878	IF(lwp) WRITE(numout,'(a, i5, a, i5, a, i4, a, i4, a, i9, a, i5, a, i5, a)') &
	879	& 'nb_cores ' , inbij,' oce + ', inbi0(ji)*inbj0(ji) - inbij &
	880	& , ' land ( ', inbi0(ji),' x ', inbj0(ji), &
	881	& ' ), nb_points ', iszi0(ji)*iszj0(ji),' ( ', iszi0(ji),' x ', iszj0(ji),' )'
	882	END DO
	883	DEALLOCATE( inbi0, inbj0, iszi0, iszj0 )
	884	RETURN
	885	ENDIF
	886
	887	DEALLOCATE( iszi0, iszj0 )
	888	inbij = inbijmax + 1 ! default: larger than possible
	889	ii = isz0+1 ! start from the end of the list (smaller subdomains)
	890	DO WHILE( inbij > knbij ) ! while the number of ocean subdomains exceed the number of procs
	891	ii = ii -1
	892	ALLOCATE( llisoce(inbi0(ii), inbj0(ii)) )
	893	CALL mpp_init_isoce( inbi0(ii), inbj0(ii), llisoce ) ! must be done by all core
	894	inbij = COUNT(llisoce)
	895	DEALLOCATE( llisoce )
	896	END DO
	897	knbi = inbi0(ii)
	898	knbj = inbj0(ii)
	899	DEALLOCATE( inbi0, inbj0 )
	900	!
	901	END SUBROUTINE mpp_init_bestpartition
	902
	903
	904	SUBROUTINE mpp_init_landprop( propland )
	905	!!----------------------------------------------------------------------
	906	!! * ROUTINE mpp_init_landprop *
	907	!!
	908	!! ** Purpose : the the proportion of land points in the surface land-sea mask
	909	!!
	910	!! ** Method : read iproc strips (of length jpiglo) of the land-sea mask
	911	!!----------------------------------------------------------------------
	912	REAL(wp), INTENT( out) :: propland ! proportion of land points (between 0 and 1)
	913	!
	914	INTEGER, DIMENSION(jpni*jpnj) :: kusedom_1d
	915	INTEGER :: inboce
	916	INTEGER :: iproc, idiv, ijsz
	917	INTEGER :: ijstr, ijend, ijcnt
	918	LOGICAL, ALLOCATABLE, DIMENSION(:,:) :: lloce
	919	!!----------------------------------------------------------------------
	920	! do nothing if there is no land-sea mask
	921	IF( numbot == -1 .and. numbdy == -1 ) THEN
	922	propland = 0.
	923	RETURN
	924	ENDIF
	925
	926	! number of processes reading the bathymetry file
	927	iproc = MINVAL( (/mppsize, jpjglo/2, 100/) ) ! read a least 2 lines, no more that 100 processes reading at the same time
	928
	929	! we want to read iproc strips of the land-sea mask. -> pick up iproc processes among mppsize processes
	930	IF( iproc == 1 ) THEN ; idiv = mppsize
	931	ELSE ; idiv = ( mppsize - 1 ) / ( iproc - 1 )
	932	ENDIF
	933	ijsz = jpjglo / iproc
	934	IF( narea <= MOD(jpjglo,iproc) ) ijsz = ijsz + 1
	935
	936	IF( MOD( narea-1, idiv ) == 0 .AND. (idiv /= 1 .OR. narea <= iproc ) ) THEN
	937	!
	938	ijstr = (narea-1)*(jpjglo/iproc) + MIN(narea-1, MOD(jpjglo,iproc)) + 1
	939	ijend = ijstr + ijsz - 1
	940	ijcnt = ijend - ijstr + 1
	941	!
	942	ALLOCATE( lloce(jpiglo, ijcnt) ) ! allocate the strip
	943	CALL mpp_init_readbot_strip( ijstr, ijcnt, lloce )
	944	inboce = COUNT(lloce)
	945	DEALLOCATE(lloce)
	946	!
	947	ELSE
	948	inboce = 0
	949	ENDIF
	950	CALL mpp_sum( 'mppini', inboce )
	951	!
	952	propland = REAL( jpiglojpjglo - inboce, wp ) / REAL( jpiglojpjglo, wp )
	953	!
	954	END SUBROUTINE mpp_init_landprop
	955
	956
	957	SUBROUTINE mpp_init_isoce( knbi, knbj, ldisoce )
	958	!!----------------------------------------------------------------------
	959	!! * ROUTINE mpp_init_nboce *
	960	!!
	961	!! ** Purpose : check for a mpi domain decomposition knbi x knbj which
	962	!! subdomains contain at least 1 ocean point
	963	!!
	964	!! ** Method : read knbj strips (of length jpiglo) of the land-sea mask
	965	!!----------------------------------------------------------------------
	966	INTEGER, INTENT(in ) :: knbi, knbj
	967	LOGICAL, DIMENSION(knbi,knbj), INTENT( out) :: ldisoce !
	968	!
	969	INTEGER, DIMENSION(knbi,knbj) :: inboce
	970	INTEGER, DIMENSION(knbi*knbj) :: inboce_1d
	971	INTEGER :: idiv, i2read, inj
	972	INTEGER :: iimax, ijmax
	973	INTEGER :: ji,jj
	974	LOGICAL, ALLOCATABLE, DIMENSION(:,:) :: lloce
	975	INTEGER, ALLOCATABLE, DIMENSION(:,:) :: iimppt, ilci
	976	INTEGER, ALLOCATABLE, DIMENSION(:,:) :: ijmppt, ilcj
	977	!!----------------------------------------------------------------------
	978	! do nothing if there is no land-sea mask
	979	IF( numbot == -1 .AND. numbdy == -1 ) THEN
	980	ldisoce(:,:) = .TRUE.
	981	RETURN
	982	ENDIF
	983
	984	! we want to read knbj strips of the land-sea mask. -> pick up knbj processes among mppsize processes
	985	IF( knbj == 1 ) THEN ; idiv = mppsize
	986	ELSE ; idiv = ( mppsize - 1 ) / ( knbj - 1 )
	987	ENDIF
	988	inboce(:,:) = 0
	989	IF( MOD( narea-1, idiv ) == 0 .AND. (idiv /= 1 .OR. narea <= knbj ) ) THEN
	990	!
	991	ALLOCATE( iimppt(knbi,knbj), ijmppt(knbi,knbj), ilci(knbi,knbj), ilcj(knbi,knbj) )
	992	CALL mpp_basic_decomposition( knbi, knbj, iimax, ijmax, iimppt, ijmppt, ilci, ilcj )
	993	!
	994	i2read = knbj / mppsize ! strip number to be read by this process
	995	IF( ( narea - 1 ) / idiv < MOD(knbj,mppsize) ) i2read = i2read + 1
	996	DO jj = 1, i2read
	997	! strip number to be read (from 1 to knbj)
	998	inj = ( narea - 1 ) * ( knbj / mppsize ) + MIN( MOD(knbj,mppsize), ( narea - 1 ) / idiv ) + jj
	999	ALLOCATE( lloce(jpiglo, ilcj(1,inj)) ) ! allocate the strip
	1000	CALL mpp_init_readbot_strip( ijmppt(1,inj), ilcj(1,inj), lloce ) ! read the strip
	1001	DO ji = 1, knbi
	1002	inboce(ji,inj) = COUNT( lloce(iimppt(ji,1):iimppt(ji,1)+ilci(ji,1)-1,:) )
	1003	END DO
	1004	DEALLOCATE(lloce)
	1005	END DO
	1006	!
	1007	DEALLOCATE(iimppt, ijmppt, ilci, ilcj)
	1008	!
	1009	ENDIF
	1010
	1011	inboce_1d = RESHAPE(inboce, (/ knbi*knbj /))
	1012	CALL mpp_sum( 'mppini', inboce_1d )
	1013	inboce = RESHAPE(inboce_1d, (/knbi, knbj/))
	1014	ldisoce = inboce /= 0
	1015	!
	1016	END SUBROUTINE mpp_init_isoce
	1017
	1018
	1019	SUBROUTINE mpp_init_readbot_strip( kjstr, kjcnt, ldoce )
	1020	!!----------------------------------------------------------------------
	1021	!! * ROUTINE mpp_init_readbot_strip *
	1022	!!
	1023	!! ** Purpose : Read relevant bathymetric information in order to
	1024	!! provide a land/sea mask used for the elimination
	1025	!! of land domains, in an mpp computation.
	1026	!!
	1027	!! ** Method : read stipe of size (jpiglo,...)
	1028	!!----------------------------------------------------------------------
	1029	INTEGER , INTENT(in ) :: kjstr !
	1030	INTEGER , INTENT(in ) :: kjcnt !
	1031	LOGICAL, DIMENSION(jpiglo,kjcnt), INTENT( out) :: ldoce !
	1032	!
	1033	INTEGER :: inumsave ! local logical unit
	1034	REAL(wp), DIMENSION(jpiglo,kjcnt) :: zbot, zbdy
	1035	!!----------------------------------------------------------------------
	1036	!
	1037	inumsave = numout ; numout = numnul ! redirect all print to /dev/null
	1038	!
	1039	IF( numbot /= -1 ) THEN
	1040	CALL iom_get( numbot, jpdom_unknown, 'bottom_level', zbot, kstart = (/1,kjstr/), kcount = (/jpiglo, kjcnt/) )
	1041	ELSE
	1042	zbot(:,:) = 1. ! put a non-null value
	1043	ENDIF
	1044
	1045	IF( numbdy /= -1 ) THEN ! Adjust with bdy_msk if it exists
	1046	CALL iom_get ( numbdy, jpdom_unknown, 'bdy_msk', zbdy, kstart = (/1,kjstr/), kcount = (/jpiglo, kjcnt/) )
	1047	zbot(:,:) = zbot(:,:) * zbdy(:,:)
	1048	ENDIF
	1049	!
	1050	ldoce = zbot > 0.
	1051	numout = inumsave
	1052	!
	1053	END SUBROUTINE mpp_init_readbot_strip
	1054
	1055
[88]	1056	SUBROUTINE mpp_init_ioipsl
	1057	!!----------------------------------------------------------------------
	1058	!! * ROUTINE mpp_init_ioipsl *
	1059	!!
	1060	!! ** Purpose :
	1061	!!
	1062	!! ** Method :
	1063	!!
	1064	!! History :
[1238]	1065	!! 9.0 ! 04-03 (G. Madec ) MPP-IOIPSL
	1066	!! " " ! 08-12 (A. Coward) addition in case of jpni*jpnj < jpnij
[88]	1067	!!----------------------------------------------------------------------
[2715]	1068	INTEGER, DIMENSION(2) :: iglo, iloc, iabsf, iabsl, ihals, ihale, idid
[88]	1069	!!----------------------------------------------------------------------
[352]	1070
[1238]	1071	! The domain is split only horizontally along i- or/and j- direction
	1072	! So we need at the most only 1D arrays with 2 elements.
	1073	! Set idompar values equivalent to the jpdom_local_noextra definition
	1074	! used in IOM. This works even if jpnij .ne. jpni*jpnj.
[88]	1075	iglo(1) = jpiglo
	1076	iglo(2) = jpjglo
	1077	iloc(1) = nlci
	1078	iloc(2) = nlcj
	1079	iabsf(1) = nimppt(narea)
	1080	iabsf(2) = njmppt(narea)
	1081	iabsl(:) = iabsf(:) + iloc(:) - 1
[1238]	1082	ihals(1) = nldi - 1
	1083	ihals(2) = nldj - 1
	1084	ihale(1) = nlci - nlei
	1085	ihale(2) = nlcj - nlej
[352]	1086	idid(1) = 1
	1087	idid(2) = 2
	1088
[88]	1089	IF(lwp) THEN
[516]	1090	WRITE(numout,*)
[352]	1091	WRITE(numout,*) 'mpp_init_ioipsl : iloc = ', iloc (1), iloc (2)
	1092	WRITE(numout,*) '~~~~~~~~~~~~~~~ iabsf = ', iabsf(1), iabsf(2)
	1093	WRITE(numout,*) ' ihals = ', ihals(1), ihals(2)
	1094	WRITE(numout,*) ' ihale = ', ihale(1), ihale(2)
[88]	1095	ENDIF
[2715]	1096	!
[352]	1097	CALL flio_dom_set ( jpnij, nproc, idid, iglo, iloc, iabsf, iabsl, ihals, ihale, 'BOX', nidom)
[2715]	1098	!
[88]	1099	END SUBROUTINE mpp_init_ioipsl
	1100
[9436]	1101
	1102	SUBROUTINE mpp_init_nfdcom
	1103	!!----------------------------------------------------------------------
	1104	!! * ROUTINE mpp_init_nfdcom *
	1105	!! ** Purpose : Setup for north fold exchanges with explicit
	1106	!! point-to-point messaging
	1107	!!
	1108	!! ** Method : Initialization of the northern neighbours lists.
	1109	!!----------------------------------------------------------------------
	1110	!! 1.0 ! 2011-10 (A. C. Coward, NOCS & J. Donners, PRACE)
	1111	!! 2.0 ! 2013-06 Setup avoiding MPI communication (I. Epicoco, S. Mocavero, CMCC)
	1112	!!----------------------------------------------------------------------
	1113	INTEGER :: sxM, dxM, sxT, dxT, jn
	1114	INTEGER :: njmppmax
	1115	!!----------------------------------------------------------------------
	1116	!
	1117	njmppmax = MAXVAL( njmppt )
	1118	!
	1119	!initializes the north-fold communication variables
	1120	isendto(:) = 0
	1121	nsndto = 0
	1122	!
	1123	IF ( njmpp == njmppmax ) THEN ! if I am a process in the north
	1124	!
	1125	!sxM is the first point (in the global domain) needed to compute the north-fold for the current process
	1126	sxM = jpiglo - nimppt(narea) - nlcit(narea) + 1
	1127	!dxM is the last point (in the global domain) needed to compute the north-fold for the current process
	1128	dxM = jpiglo - nimppt(narea) + 2
	1129	!
	1130	! loop over the other north-fold processes to find the processes
	1131	! managing the points belonging to the sxT-dxT range
	1132	!
	1133	DO jn = 1, jpni
	1134	!
	1135	sxT = nfiimpp(jn, jpnj) ! sxT = 1st point (in the global domain) of the jn process
	1136	dxT = nfiimpp(jn, jpnj) + nfilcit(jn, jpnj) - 1 ! dxT = last point (in the global domain) of the jn process
	1137	!
	1138	IF ( sxT < sxM .AND. sxM < dxT ) THEN
	1139	nsndto = nsndto + 1
	1140	isendto(nsndto) = jn
	1141	ELSEIF( sxM <= sxT .AND. dxM >= dxT ) THEN
	1142	nsndto = nsndto + 1
	1143	isendto(nsndto) = jn
	1144	ELSEIF( dxM < dxT .AND. sxT < dxM ) THEN
	1145	nsndto = nsndto + 1
	1146	isendto(nsndto) = jn
	1147	ENDIF
	1148	!
	1149	END DO
	1150	nfsloop = 1
	1151	nfeloop = nlci
	1152	DO jn = 2,jpni-1
	1153	IF( nfipproc(jn,jpnj) == (narea - 1) ) THEN
	1154	IF( nfipproc(jn-1,jpnj) == -1 ) nfsloop = nldi
	1155	IF( nfipproc(jn+1,jpnj) == -1 ) nfeloop = nlei
	1156	ENDIF
	1157	END DO
	1158	!
	1159	ENDIF
	1160	l_north_nogather = .TRUE.
	1161	!
	1162	END SUBROUTINE mpp_init_nfdcom
	1163
[10330]	1164
[3]	1165	#endif
[88]	1166
[3]	1167	!!======================================================================
	1168	END MODULE mppini

Note: See TracBrowser for help on using the repository browser.

Download in other formats: