Context Navigation

← Previous Revision
Next Revision →
Normal
Revision Log

mppini.F90

Last change on this file was 15302, checked in by smasson, 3 years ago
trunk: minor improvment in pt2pt (east-west comm only send/recv inner values)
Property svn:keywords set to `Id`
File size: 71.6 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
[13286]	10	!! 3.4 ! 2011-10 (A. C. Coward, NOCS & J. Donners, PRACE) add init_nfdcom
[14072]	11	!! 3. ! 2013-06 (I. Epicoco, S. Mocavero, CMCC) 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	!!----------------------------------------------------------------------
[13286]	17	!! mpp_init : Lay out the global domain over processors with/without land processor elimination
[14072]	18	!! init_ioipsl: IOIPSL initialization in mpp
[13286]	19	!! init_nfdcom: Setup for north fold exchanges with explicit point-to-point messaging
[14072]	20	!! init_doloop: set the starting/ending indices of DO-loop used in do_loop_substitute
[3]	21	!!----------------------------------------------------------------------
[9019]	22	USE dom_oce ! ocean space and time domain
[14072]	23	USE bdy_oce ! open BounDarY
[9019]	24	!
[15267]	25	USE lbcnfd ! Setup of north fold exchanges
[9019]	26	USE lib_mpp ! distribued memory computing library
[14072]	27	USE iom ! nemo I/O library
[9019]	28	USE ioipsl ! I/O IPSL library
	29	USE in_out_manager ! I/O Manager
[3]	30
	31	IMPLICIT NONE
	32	PRIVATE
	33
[13286]	34	PUBLIC mpp_init ! called by nemogcm.F90
	35	PUBLIC mpp_getnum ! called by prtctl
	36	PUBLIC mpp_basesplit ! called by prtctl
	37	PUBLIC mpp_is_ocean ! called by prtctl
[14072]	38
[13286]	39	INTEGER :: numbot = -1 ! 'bottom_level' local logical unit
	40	INTEGER :: numbdy = -1 ! 'bdy_msk' local logical unit
[14072]	41
[3]	42	!!----------------------------------------------------------------------
[9598]	43	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[14072]	44	!! $Id$
[10068]	45	!! Software governed by the CeCILL license (see ./LICENSE)
[3]	46	!!----------------------------------------------------------------------
	47	CONTAINS
	48
[14229]	49	#if defined key_mpi_off
[3]	50	!!----------------------------------------------------------------------
	51	!! Default option : shared memory computing
	52	!!----------------------------------------------------------------------
	53
	54	SUBROUTINE mpp_init
	55	!!----------------------------------------------------------------------
	56	!! * ROUTINE mpp_init *
	57	!!
	58	!! ** Purpose : Lay out the global domain over processors.
	59	!!
	60	!! ** Method : Shared memory computing, set the local processor
[9019]	61	!! variables to the value of the global domain
[3]	62	!!----------------------------------------------------------------------
[9019]	63	!
[15023]	64	nn_comm = 1
	65	nn_hls = 1
	66	jpiglo = Ni0glo + 2 * nn_hls
	67	jpjglo = Nj0glo + 2 * nn_hls
	68	jpimax = jpiglo
	69	jpjmax = jpjglo
	70	jpi = jpiglo
	71	jpj = jpjglo
	72	jpk = MAX( 2, jpkglo )
	73	jpij = jpi*jpj
	74	jpni = 1
	75	jpnj = 1
	76	jpnij = jpni*jpnj
	77	nimpp = 1
	78	njmpp = 1
	79	nidom = FLIO_DOM_NONE
[9019]	80	!
[15023]	81	mpiSnei(:,:) = -1
	82	mpiRnei(:,:) = -1
	83	l_SelfPerio(1:2) = l_Iperio ! west, east periodicity by itself
	84	l_SelfPerio(3:4) = l_Jperio ! south, north periodicity by itself
	85	l_SelfPerio(5:8) = l_Iperio .AND. l_Jperio ! corners bi-periodicity by itself
	86	l_IdoNFold = l_NFold ! is this process doing North fold?
	87	!
[14072]	88	CALL init_doloop ! set start/end indices or do-loop depending on the halo width value (nn_hls)
[13438]	89	!
[3]	90	IF(lwp) THEN
	91	WRITE(numout,*)
[9019]	92	WRITE(numout,*) 'mpp_init : NO massively parallel processing'
	93	WRITE(numout,*) '~~~~~~~~ '
[14072]	94	WRITE(numout,*) ' l_Iperio = ', l_Iperio, ' l_Jperio = ', l_Jperio
[14433]	95	WRITE(numout,*) ' njmpp = ', njmpp
[3]	96	ENDIF
[9019]	97	!
[13216]	98	#if defined key_agrif
	99	call agrif_nemo_init()
	100	#endif
[3]	101	END SUBROUTINE mpp_init
	102
	103	#else
	104	!!----------------------------------------------------------------------
[14229]	105	!! MPI massively parallel processing
[3]	106	!!----------------------------------------------------------------------
	107
[9449]	108
[3]	109	SUBROUTINE mpp_init
	110	!!----------------------------------------------------------------------
	111	!! * ROUTINE mpp_init *
[14072]	112	!!
[3]	113	!! ** Purpose : Lay out the global domain over processors.
[9019]	114	!! If land processors are to be eliminated, this program requires the
	115	!! presence of the domain configuration file. Land processors elimination
	116	!! is performed if jpni x jpnj /= jpnij. In this case, using the MPP_PREP
	117	!! preprocessing tool, help for defining the best cutting out.
[3]	118	!!
	119	!! ** Method : Global domain is distributed in smaller local domains.
	120	!! Periodic condition is a function of the local domain position
[14433]	121	!! (global boundary or neighbouring domain) and of the global periodic
[3]	122	!!
[9019]	123	!! ** Action : - set domain parameters
[14072]	124	!! nimpp : longitudinal index
[3]	125	!! njmpp : latitudinal index
	126	!! narea : number for local area
[14433]	127	!! mpinei : number of neighboring domains (starting at 0, -1 if no neighbourg)
[3]	128	!!----------------------------------------------------------------------
[14433]	129	INTEGER :: ji, jj, jn, jp, jh
	130	INTEGER :: ii, ij, ii2, ij2
	131	INTEGER :: inijmin ! number of oce subdomains
	132	INTEGER :: inum, inum0
	133	INTEGER :: ifreq, il1, imil, il2, ijm1
	134	INTEGER :: ierr, ios
	135	INTEGER :: inbi, inbj, iimax, ijmax, icnt1, icnt2
	136	INTEGER, DIMENSION(16*n_hlsmax) :: ichanged
	137	INTEGER, ALLOCATABLE, DIMENSION(: ) :: iin, ijn
	138	INTEGER, ALLOCATABLE, DIMENSION(:,: ) :: iimppt, ijpi, ipproc
	139	INTEGER, ALLOCATABLE, DIMENSION(:,: ) :: ijmppt, ijpj
	140	INTEGER, ALLOCATABLE, DIMENSION(:,: ) :: impi
	141	INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: inei
[10615]	142	LOGICAL :: llbest, llauto
[10570]	143	LOGICAL :: llwrtlay
[14433]	144	LOGICAL :: llmpi_Iperio, llmpi_Jperio, llmpiNFold
[11536]	145	LOGICAL :: ln_listonly
[14433]	146	LOGICAL, ALLOCATABLE, DIMENSION(:,: ) :: llisOce ! is not land-domain only?
	147	LOGICAL, ALLOCATABLE, DIMENSION(:,:,:) :: llnei ! are neighbourgs existing?
[10425]	148	NAMELIST/nambdy/ ln_bdy, nb_bdy, ln_coords_file, cn_coords_file, &
	149	& ln_mask_file, cn_mask_file, cn_dyn2d, nn_dyn2d_dta, &
[14072]	150	& cn_dyn3d, nn_dyn3d_dta, cn_tra, nn_tra_dta, &
[10425]	151	& ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp, rn_time_dmp_out, &
	152	& cn_ice, nn_ice_dta, &
[11536]	153	& ln_vol, nn_volctl, nn_rimwidth
[14433]	154	NAMELIST/nammpp/ jpni, jpnj, nn_hls, ln_nnogather, ln_listonly, nn_comm
[3]	155	!!----------------------------------------------------------------------
[11536]	156	!
[12377]	157	llwrtlay = lwm .OR. sn_cfctl%l_layout
[11536]	158	!
	159	! 0. read namelists parameters
	160	! -----------------------------------
	161	!
	162	READ ( numnam_ref, nammpp, IOSTAT = ios, ERR = 901 )
	163	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nammpp in reference namelist' )
	164	READ ( numnam_cfg, nammpp, IOSTAT = ios, ERR = 902 )
[14072]	165	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'nammpp in configuration namelist' )
[11536]	166	!
[13286]	167	nn_hls = MAX(1, nn_hls) ! nn_hls must be > 0
[11536]	168	IF(lwp) THEN
	169	WRITE(numout,*) ' Namelist nammpp'
[15296]	170	IF( jpni < 1 .OR. jpnj < 1 ) THEN
[11536]	171	WRITE(numout,*) ' jpni and jpnj will be calculated automatically'
	172	ELSE
	173	WRITE(numout,*) ' processor grid extent in i jpni = ', jpni
	174	WRITE(numout,*) ' processor grid extent in j jpnj = ', jpnj
	175	ENDIF
	176	WRITE(numout,*) ' avoid use of mpi_allgather at the north fold ln_nnogather = ', ln_nnogather
[13286]	177	WRITE(numout,*) ' halo width (applies to both rows and columns) nn_hls = ', nn_hls
[11536]	178	ENDIF
	179	!
	180	IF(lwm) WRITE( numond, nammpp )
[13286]	181	!
	182	jpiglo = Ni0glo + 2 * nn_hls
	183	jpjglo = Nj0glo + 2 * nn_hls
	184	!
[10425]	185	! do we need to take into account bdy_msk?
	186	READ ( numnam_ref, nambdy, IOSTAT = ios, ERR = 903)
[11536]	187	903 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy in reference namelist (mppini)' )
[10425]	188	READ ( numnam_cfg, nambdy, IOSTAT = ios, ERR = 904 )
[11536]	189	904 IF( ios > 0 ) CALL ctl_nam ( ios , 'nambdy in configuration namelist (mppini)' )
[10425]	190	!
	191	IF( ln_read_cfg ) CALL iom_open( cn_domcfg, numbot )
	192	IF( ln_bdy .AND. ln_mask_file ) CALL iom_open( cn_mask_file, numbdy )
	193	!
[13286]	194	IF( ln_listonly ) CALL bestpartition( MAX(mppsize,jpni*jpnj), ldlist = .TRUE. ) ! must be done by all core
[11536]	195	!
[10425]	196	! 1. Dimension arrays for subdomains
	197	! -----------------------------------
	198	!
[14433]	199	! If dimensions of MPI processes grid weren't specified in the namelist file
[9436]	200	! then we calculate them here now that we have our communicator size
[10615]	201	IF(lwp) THEN
[14053]	202	WRITE(numout,*)
[10615]	203	WRITE(numout,*) 'mpp_init:'
	204	WRITE(numout,*) '~~~~~~~~ '
	205	ENDIF
[10425]	206	IF( jpni < 1 .OR. jpnj < 1 ) THEN
[13286]	207	CALL bestpartition( mppsize, jpni, jpnj ) ! best mpi decomposition for mppsize mpi processes
[10615]	208	llauto = .TRUE.
[10425]	209	llbest = .TRUE.
	210	ELSE
[10615]	211	llauto = .FALSE.
[13286]	212	CALL bestpartition( mppsize, inbi, inbj, icnt2 ) ! best mpi decomposition for mppsize mpi processes
[10615]	213	! largest subdomain size for mpi decoposition jpni*jpnj given in the namelist
[13286]	214	CALL mpp_basesplit( jpiglo, jpjglo, nn_hls, jpni, jpnj, jpimax, jpjmax )
	215	! largest subdomain size for mpi decoposition inbi*inbj given by bestpartition
	216	CALL mpp_basesplit( jpiglo, jpjglo, nn_hls, inbi, inbj, iimax, ijmax )
[10615]	217	icnt1 = jpni*jpnj - mppsize ! number of land subdomains that should be removed to use mppsize mpi processes
	218	IF(lwp) THEN
	219	WRITE(numout,9000) ' The chosen domain decomposition ', jpni, ' x ', jpnj, ' with ', icnt1, ' land subdomains'
	220	WRITE(numout,9002) ' - uses a total of ', mppsize,' mpi process'
	221	WRITE(numout,9000) ' - has mpi subdomains with a maximum size of (jpi = ', jpimax, ', jpj = ', jpjmax, &
	222	& ', jpijpj = ', jpimaxjpjmax, ')'
	223	WRITE(numout,9000) ' The best domain decompostion ', inbi, ' x ', inbj, ' with ', icnt2, ' land subdomains'
	224	WRITE(numout,9002) ' - uses a total of ', inbi*inbj-icnt2,' mpi process'
	225	WRITE(numout,9000) ' - has mpi subdomains with a maximum size of (jpi = ', iimax, ', jpj = ', ijmax, &
	226	& ', jpijpj = ', iimax ijmax, ')'
	227	ENDIF
	228	IF( iimaxijmax < jpimaxjpjmax ) THEN ! chosen subdomain size is larger that the best subdomain size
[10425]	229	llbest = .FALSE.
[10615]	230	IF ( inbi*inbj-icnt2 < mppsize ) THEN
	231	WRITE(ctmp1,*) ' ==> You could therefore have smaller mpi subdomains with less mpi processes'
	232	ELSE
	233	WRITE(ctmp1,*) ' ==> You could therefore have smaller mpi subdomains with the same number of mpi processes'
	234	ENDIF
	235	CALL ctl_warn( ' ', ctmp1, ' ', ' --- YOU ARE WASTING CPU... ---', ' ' )
	236	ELSE IF ( iimaxijmax == jpimaxjpjmax .AND. (inbi*inbj-icnt2) < mppsize) THEN
	237	llbest = .FALSE.
	238	WRITE(ctmp1,*) ' ==> You could therefore have the same mpi subdomains size with less mpi processes'
	239	CALL ctl_warn( ' ', ctmp1, ' ', ' --- YOU ARE WASTING CPU... ---', ' ' )
[10425]	240	ELSE
	241	llbest = .TRUE.
	242	ENDIF
	243	ENDIF
[14072]	244
[10425]	245	! look for land mpi subdomains...
[14433]	246	ALLOCATE( llisOce(jpni,jpnj) )
	247	CALL mpp_is_ocean( llisOce )
	248	inijmin = COUNT( llisOce ) ! number of oce subdomains
[10425]	249
[10615]	250	IF( mppsize < inijmin ) THEN ! too many oce subdomains: can happen only if jpni and jpnj are prescribed...
[10425]	251	WRITE(ctmp1,9001) ' With this specified domain decomposition: jpni = ', jpni, ' jpnj = ', jpnj
	252	WRITE(ctmp2,9002) ' we can eliminate only ', jpni*jpnj - inijmin, ' land mpi subdomains therefore '
	253	WRITE(ctmp3,9001) ' the number of ocean mpi subdomains (', inijmin,') exceed the number of MPI processes:', mppsize
	254	WRITE(ctmp4,*) ' ==>>> There is the list of best domain decompositions you should use: '
[10615]	255	CALL ctl_stop( ctmp1, ctmp2, ctmp3, ' ', ctmp4, ' ' )
[13286]	256	CALL bestpartition( mppsize, ldlist = .TRUE. ) ! must be done by all core
[10425]	257	ENDIF
	258
[10615]	259	IF( mppsize > jpni*jpnj ) THEN ! not enough mpi subdomains for the total number of mpi processes
	260	IF(lwp) THEN
	261	WRITE(numout,9003) ' The number of mpi processes: ', mppsize
	262	WRITE(numout,9003) ' exceeds the maximum number of subdomains (ocean+land) = ', jpni*jpnj
	263	WRITE(numout,9001) ' defined by the following domain decomposition: jpni = ', jpni, ' jpnj = ', jpnj
	264	WRITE(numout, *) ' You should: '
	265	IF( llauto ) THEN
	266	WRITE(numout,*) ' - either prescribe your domain decomposition with the namelist variables'
	267	WRITE(numout,*) ' jpni and jpnj to match the number of mpi process you want to use, '
	268	WRITE(numout,*) ' even IF it not the best choice...'
	269	WRITE(numout,*) ' - or keep the automatic and optimal domain decomposition by picking up one'
	270	WRITE(numout,*) ' of the number of mpi process proposed in the list bellow'
	271	ELSE
	272	WRITE(numout,*) ' - either properly prescribe your domain decomposition with jpni and jpnj'
	273	WRITE(numout,*) ' in order to be consistent with the number of mpi process you want to use'
	274	WRITE(numout,*) ' even IF it not the best choice...'
	275	WRITE(numout,*) ' - or use the automatic and optimal domain decomposition and pick up one of'
	276	WRITE(numout,*) ' the domain decomposition proposed in the list bellow'
	277	ENDIF
	278	WRITE(numout,*)
	279	ENDIF
[13286]	280	CALL bestpartition( mppsize, ldlist = .TRUE. ) ! must be done by all core
[10425]	281	ENDIF
	282
	283	jpnij = mppsize ! force jpnij definition <-- remove as much land subdomains as needed to reach this condition
	284	IF( mppsize > inijmin ) THEN
	285	WRITE(ctmp1,9003) ' The number of mpi processes: ', mppsize
	286	WRITE(ctmp2,9003) ' exceeds the maximum number of ocean subdomains = ', inijmin
	287	WRITE(ctmp3,9002) ' we suppressed ', jpni*jpnj - mppsize, ' land subdomains '
	288	WRITE(ctmp4,9002) ' BUT we had to keep ', mppsize - inijmin, ' land subdomains that are useless...'
[10615]	289	CALL ctl_warn( ctmp1, ctmp2, ctmp3, ctmp4, ' ', ' --- YOU ARE WASTING CPU... ---', ' ' )
[10425]	290	ELSE ! mppsize = inijmin
	291	IF(lwp) THEN
[10615]	292	IF(llbest) WRITE(numout,*) ' ==> you use the best mpi decomposition'
	293	WRITE(numout,*)
[10425]	294	WRITE(numout,9003) ' Number of mpi processes: ', mppsize
	295	WRITE(numout,9003) ' Number of ocean subdomains = ', inijmin
	296	WRITE(numout,9003) ' Number of suppressed land subdomains = ', jpni*jpnj - inijmin
	297	WRITE(numout,*)
	298	ENDIF
	299	ENDIF
	300	9000 FORMAT (a, i4, a, i4, a, i7, a)
	301	9001 FORMAT (a, i4, a, i4)
	302	9002 FORMAT (a, i4, a)
	303	9003 FORMAT (a, i5)
	304
[14433]	305	ALLOCATE( nfimpp(jpni), nfproc(jpni), nfjpi(jpni), &
	306	& iin(jpnij), ijn(jpnij), &
	307	& iimppt(jpni,jpnj), ijmppt(jpni,jpnj), ijpi(jpni,jpnj), ijpj(jpni,jpnj), ipproc(jpni,jpnj), &
	308	& inei(8,jpni,jpnj), llnei(8,jpni,jpnj), &
	309	& impi(8,jpnij), &
	310	& STAT=ierr )
[10425]	311	CALL mpp_sum( 'mppini', ierr )
[9436]	312	IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'mpp_init: unable to allocate standard ocean arrays' )
[14072]	313
[9508]	314	#if defined key_agrif
[13216]	315	CALL agrif_nemo_init()
[9508]	316	#endif
[9436]	317	!
[10425]	318	! 2. Index arrays for subdomains
[3]	319	! -----------------------------------
[9019]	320	!
[13286]	321	CALL mpp_basesplit( jpiglo, jpjglo, nn_hls, jpni, jpnj, jpimax, jpjmax, iimppt, ijmppt, ijpi, ijpj )
[14433]	322	CALL mpp_getnum( llisOce, ipproc, iin, ijn )
[9019]	323	!
[14433]	324	ii = iin(narea)
	325	ij = ijn(narea)
	326	jpi = ijpi(ii,ij)
	327	jpj = ijpj(ii,ij)
	328	jpk = MAX( 2, jpkglo )
	329	jpij = jpi*jpj
	330	nimpp = iimppt(ii,ij)
	331	njmpp = ijmppt(ii,ij)
[13286]	332	!
[15267]	333	CALL init_doloop ! set start/end indices of do-loop, depending on the halo width value (nn_hls)
	334	CALL init_locglo ! define now functions needed to convert indices from/to global to/from local domains
[14433]	335	!
[9019]	336	IF(lwp) THEN
	337	WRITE(numout,*)
[10425]	338	WRITE(numout,*) 'MPI Message Passing MPI - domain lay out over processors'
	339	WRITE(numout,*)
[9019]	340	WRITE(numout,*) ' defines mpp subdomains'
[14072]	341	WRITE(numout,*) ' jpni = ', jpni
[10425]	342	WRITE(numout,*) ' jpnj = ', jpnj
[13286]	343	WRITE(numout,*) ' jpnij = ', jpnij
[14433]	344	WRITE(numout,*) ' nimpp = ', nimpp
	345	WRITE(numout,*) ' njmpp = ', njmpp
[9019]	346	WRITE(numout,*)
[13286]	347	WRITE(numout,*) ' sum ijpi(i,1) = ', sum(ijpi(:,1)), ' jpiglo = ', jpiglo
[14433]	348	WRITE(numout,*) ' sum ijpj(1,j) = ', SUM(ijpj(1,:)), ' jpjglo = ', jpjglo
	349
	350	! Subdomain grid print
[10425]	351	ifreq = 4
	352	il1 = 1
	353	DO jn = 1, (jpni-1)/ifreq+1
	354	il2 = MIN(jpni,il1+ifreq-1)
	355	WRITE(numout,*)
	356	WRITE(numout,9400) ('***',ji=il1,il2-1)
	357	DO jj = jpnj, 1, -1
	358	WRITE(numout,9403) (' ',ji=il1,il2-1)
[13286]	359	WRITE(numout,9402) jj, (ijpi(ji,jj),ijpj(ji,jj),ji=il1,il2)
[10425]	360	WRITE(numout,9404) (ipproc(ji,jj),ji=il1,il2)
	361	WRITE(numout,9403) (' ',ji=il1,il2-1)
	362	WRITE(numout,9400) ('***',ji=il1,il2-1)
	363	END DO
	364	WRITE(numout,9401) (ji,ji=il1,il2)
	365	il1 = il1+ifreq
	366	END DO
	367	9400 FORMAT(' *' ,20('***********',a3) )
	368	9403 FORMAT(' * ',20(' * ',a3) )
	369	9401 FORMAT(' ' ,20(' ',i3,' ') )
	370	9402 FORMAT(' ',i3,' * ',20(i3,' x',i3,' * ') )
[11640]	371	9404 FORMAT(' * ' ,20(' ' ,i4,' * ') )
[10425]	372	ENDIF
[14433]	373	!
	374	! Store informations for the north pole folding communications
	375	nfproc(:) = ipproc(:,jpnj)
	376	nfimpp(:) = iimppt(:,jpnj)
	377	nfjpi (:) = ijpi(:,jpnj)
	378	!
	379	! 3. Define Western, Eastern, Southern and Northern neighbors + corners in the subdomain grid reference
	380	! ------------------------------------------------------------------------------------------------------
	381	!
	382	! note that North fold is has specific treatment for its MPI communications.
	383	! This must not be treated as a "usual" communication with a northern neighbor.
	384	! -> North fold processes have no Northern neighbor in the definition done bellow
	385	!
	386	llmpi_Iperio = jpni > 1 .AND. l_Iperio ! do i-periodicity with an MPI communication?
	387	llmpi_Jperio = jpnj > 1 .AND. l_Jperio ! do j-periodicity with an MPI communication?
	388	!
	389	l_SelfPerio(1:2) = l_Iperio .AND. jpni == 1 ! west, east periodicity by itself
	390	l_SelfPerio(3:4) = l_Jperio .AND. jpnj == 1 ! south, north periodicity by itself
	391	l_SelfPerio(5:8) = l_SelfPerio(jpwe) .AND. l_SelfPerio(jpso) ! corners bi-periodicity by itself
	392	!
	393	! define neighbors mapping (1/2): default definition: ignore if neighbours are land-only subdomains or not
	394	DO jj = 1, jpnj
	395	DO ji = 1, jpni
	396	!
	397	IF ( llisOce(ji,jj) ) THEN ! this subdomain has some ocean: it has neighbours
	398	!
	399	inum0 = ji - 1 + ( jj - 1 ) * jpni ! index in the subdomains grid. start at 0
	400	!
	401	! Is there a neighbor?
	402	llnei(jpwe,ji,jj) = ji > 1 .OR. llmpi_Iperio ! West nei exists if not the first column or llmpi_Iperio
	403	llnei(jpea,ji,jj) = ji < jpni .OR. llmpi_Iperio ! East nei exists if not the last column or llmpi_Iperio
	404	llnei(jpso,ji,jj) = jj > 1 .OR. llmpi_Jperio ! South nei exists if not the first line or llmpi_Jperio
	405	llnei(jpno,ji,jj) = jj < jpnj .OR. llmpi_Jperio ! North nei exists if not the last line or llmpi_Jperio
	406	llnei(jpsw,ji,jj) = llnei(jpwe,ji,jj) .AND. llnei(jpso,ji,jj) ! So-We nei exists if both South and West nei exist
	407	llnei(jpse,ji,jj) = llnei(jpea,ji,jj) .AND. llnei(jpso,ji,jj) ! So-Ea nei exists if both South and East nei exist
	408	llnei(jpnw,ji,jj) = llnei(jpwe,ji,jj) .AND. llnei(jpno,ji,jj) ! No-We nei exists if both North and West nei exist
	409	llnei(jpne,ji,jj) = llnei(jpea,ji,jj) .AND. llnei(jpno,ji,jj) ! No-Ea nei exists if both North and East nei exist
	410	!
	411	! Which index (starting at 0) have neighbors in the subdomains grid?
	412	IF( llnei(jpwe,ji,jj) ) inei(jpwe,ji,jj) = inum0 - 1 + jpni * COUNT( (/ ji == 1 /) )
	413	IF( llnei(jpea,ji,jj) ) inei(jpea,ji,jj) = inum0 + 1 - jpni * COUNT( (/ ji == jpni /) )
	414	IF( llnei(jpso,ji,jj) ) inei(jpso,ji,jj) = inum0 - jpni + jpni * jpnj * COUNT( (/ jj == 1 /) )
	415	IF( llnei(jpno,ji,jj) ) inei(jpno,ji,jj) = inum0 + jpni - jpni * jpnj * COUNT( (/ jj == jpnj /) )
	416	IF( llnei(jpsw,ji,jj) ) inei(jpsw,ji,jj) = inei(jpso,ji,jj) - 1 + jpni * COUNT( (/ ji == 1 /) )
	417	IF( llnei(jpse,ji,jj) ) inei(jpse,ji,jj) = inei(jpso,ji,jj) + 1 - jpni * COUNT( (/ ji == jpni /) )
	418	IF( llnei(jpnw,ji,jj) ) inei(jpnw,ji,jj) = inei(jpno,ji,jj) - 1 + jpni * COUNT( (/ ji == 1 /) )
	419	IF( llnei(jpne,ji,jj) ) inei(jpne,ji,jj) = inei(jpno,ji,jj) + 1 - jpni * COUNT( (/ ji == jpni /) )
	420	!
	421	ELSE ! land-only domain has no neighbour
	422	llnei(:,ji,jj) = .FALSE.
	423	ENDIF
	424	!
	425	END DO
[9019]	426	END DO
	427	!
[14433]	428	! define neighbors mapping (2/2): check if neighbours are not land-only subdomains
	429	DO jj = 1, jpnj
	430	DO ji = 1, jpni
	431	DO jn = 1, 8
	432	IF( llnei(jn,ji,jj) ) THEN ! if a neighbour is existing -> this should not be a land-only domain
	433	ii = 1 + MOD( inei(jn,ji,jj) , jpni )
	434	ij = 1 + inei(jn,ji,jj) / jpni
	435	llnei(jn,ji,jj) = llisOce( ii, ij )
	436	ENDIF
	437	END DO
	438	END DO
	439	END DO
[13286]	440	!
[14433]	441	! update index of the neighbours in the subdomains grid
	442	WHERE( .NOT. llnei ) inei = -1
[13286]	443	!
[9019]	444	! Save processor layout in ascii file
[10570]	445	IF (llwrtlay) THEN
[7646]	446	CALL ctl_opn( inum, 'layout.dat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE., narea )
[14433]	447	WRITE(inum,'(a)') ' jpnij jpimax jpjmax jpk jpiglo jpjglo ( local: narea jpi jpj )'
	448	WRITE(inum,'(6i7,a,3i7,a)') jpnij,jpimax,jpjmax,jpk,jpiglo,jpjglo,' ( local: ',narea,jpi,jpj,' )'
	449	WRITE(inum,*)
	450	WRITE(inum, *) '------------------------------------'
	451	WRITE(inum,'(a,i2)') ' Mapping of the default neighnourgs '
	452	WRITE(inum, *) '------------------------------------'
	453	WRITE(inum,*)
	454	WRITE(inum,'(a)') ' rank ii ij jpi jpj nimpp njmpp mpiwe mpiea mpiso mpino mpisw mpise mpinw mpine'
	455	DO jp = 1, jpnij
	456	ii = iin(jp)
	457	ij = ijn(jp)
	458	WRITE(inum,'(15i6)') jp-1, ii, ij, ijpi(ii,ij), ijpj(ii,ij), iimppt(ii,ij), ijmppt(ii,ij), inei(:,ii,ij)
	459	END DO
	460	ENDIF
[9019]	461
[14433]	462	!
	463	! 4. Define Western, Eastern, Southern and Northern neighbors + corners for each mpi process
	464	! ------------------------------------------------------------------------------------------
	465	!
	466	! rewrite information from "subdomain grid" to mpi process list
	467	! Warning, for example:
	468	! position of the northern neighbor in the "subdomain grid"
	469	! position of the northern neighbor in the "mpi process list"
	470
	471	! default definition: no neighbors
	472	impi(:,:) = -1 ! (starting at 0, -1 if no neighbourg)
	473
	474	DO jp = 1, jpnij
	475	ii = iin(jp)
	476	ij = ijn(jp)
	477	DO jn = 1, 8
	478	IF( llnei(jn,ii,ij) ) THEN ! must be tested as some land-domain can be kept to fit mppsize
	479	ii2 = 1 + MOD( inei(jn,ii,ij) , jpni )
	480	ij2 = 1 + inei(jn,ii,ij) / jpni
	481	impi(jn,jp) = ipproc( ii2, ij2 )
	482	ENDIF
[7646]	483	END DO
[14433]	484	END DO
[3]	485
[9019]	486	!
[14433]	487	! 4. keep information for the local process
	488	! -----------------------------------------
	489	!
	490	! set default neighbours
	491	mpinei(:) = impi(:,narea)
	492	DO jh = 1, n_hlsmax
	493	mpiSnei(jh,:) = impi(:,narea) ! default definition
	494	mpiRnei(jh,:) = impi(:,narea)
	495	END DO
	496	!
[3]	497	IF(lwp) THEN
[6412]	498	WRITE(numout,*)
[9169]	499	WRITE(numout,*) ' resulting internal parameters : '
[14433]	500	WRITE(numout,*) ' narea = ', narea
	501	WRITE(numout,*) ' mpi nei west = ', mpinei(jpwe) , ' mpi nei east = ', mpinei(jpea)
	502	WRITE(numout,*) ' mpi nei south = ', mpinei(jpso) , ' mpi nei north = ', mpinei(jpno)
	503	WRITE(numout,*) ' mpi nei so-we = ', mpinei(jpsw) , ' mpi nei so-ea = ', mpinei(jpse)
	504	WRITE(numout,*) ' mpi nei no-we = ', mpinei(jpnw) , ' mpi nei no-ea = ', mpinei(jpne)
[3]	505	ENDIF
[14433]	506	!
	507	CALL mpp_ini_nc(nn_hls) ! Initialize communicator for neighbourhood collective communications
	508	DO jh = 1, n_hlsmax
	509	mpi_nc_com4(jh) = mpi_nc_com4(nn_hls) ! default definition
	510	mpi_nc_com8(jh) = mpi_nc_com8(nn_hls)
	511	END DO
[15267]	512	! ! Exclude exchanges which contain only land points
[13982]	513	!
[15267]	514	IF( jpnij > 1 ) CALL init_excl_landpt
[14072]	515	!
[15267]	516	! ! Prepare mpp north fold
	517	!
	518	llmpiNFold = jpni > 1 .AND. l_NFold ! is the North fold done with an MPI communication?
	519	l_IdoNFold = ijn(narea) == jpnj .AND. l_NFold ! is this process doing North fold?
	520	!
	521	IF( llmpiNFold ) CALL init_nfdcom( llwrtlay, inum ) ! init northfold communication, must be done after init_excl_landpt
	522	!
	523	! ! Save processor layout changes in ascii file
	524	!
[14433]	525	DO jh = 1, n_hlsmax ! different halo size
	526	DO ji = 1, 8
	527	ichanged(16*(jh-1) +ji) = COUNT( mpinei(ji:ji) /= mpiSnei(jh,ji:ji) )
	528	ichanged(16*(jh-1)+8+ji) = COUNT( mpinei(ji:ji) /= mpiRnei(jh,ji:ji) )
	529	END DO
	530	END DO
	531	CALL mpp_sum( "mpp_init", ichanged ) ! must be called by all processes
	532	IF (llwrtlay) THEN
	533	WRITE(inum,*)
	534	WRITE(inum, *) '----------------------------------------------------------------------'
	535	WRITE(inum,'(a,i2)') ' Mapping of the neighnourgs once excluding comm with only land points '
	536	WRITE(inum, *) '----------------------------------------------------------------------'
	537	DO jh = 1, n_hlsmax ! different halo size
	538	WRITE(inum,*)
	539	WRITE(inum,'(a,i2)') 'halo size: ', jh
	540	WRITE(inum, *) '---------'
	541	WRITE(inum,'(a)') ' rank ii ij mpiwe mpiea mpiso mpino mpisw mpise mpinw mpine'
	542	WRITE(inum, '(11i6,a)') narea-1, iin(narea), ijn(narea), mpinei(:), ' <- Org'
	543	WRITE(inum,'(18x,8i6,a,i1,a)') mpiSnei(jh,:), ' <- Send ', COUNT( mpinei(:) /= mpiSnei(jh,:) ), ' modif'
	544	WRITE(inum,'(18x,8i6,a,i1,a)') mpiRnei(jh,:), ' <- Recv ', COUNT( mpinei(:) /= mpiRnei(jh,:) ), ' modif'
	545	WRITE(inum,*) ' total changes among all mpi tasks:'
	546	WRITE(inum,*) ' mpiwe mpiea mpiso mpino mpisw mpise mpinw mpine'
	547	WRITE(inum,'(a,8i6)') ' Send: ', ichanged(jh16-15:jh16-8)
	548	WRITE(inum,'(a,8i6)') ' Recv: ', ichanged(jh16 -7:jh16 )
	549	END DO
[9917]	550	ENDIF
[9436]	551	!
[14433]	552	CALL init_ioipsl ! Prepare NetCDF output file (if necessary)
	553	!
[14072]	554	IF (llwrtlay) CLOSE(inum)
[9917]	555	!
[14433]	556	DEALLOCATE(iin, ijn, iimppt, ijmppt, ijpi, ijpj, ipproc, inei, llnei, impi, llisOce)
[9444]	557	!
[9019]	558	END SUBROUTINE mpp_init
[3]	559
[13438]	560	#endif
[3]	561
[13286]	562	SUBROUTINE mpp_basesplit( kiglo, kjglo, khls, knbi, knbj, kimax, kjmax, kimppt, kjmppt, klci, klcj)
[9019]	563	!!----------------------------------------------------------------------
[13286]	564	!! * ROUTINE mpp_basesplit *
[14072]	565	!!
[10425]	566	!! ** Purpose : Lay out the global domain over processors.
[9019]	567	!!
[10425]	568	!! ** Method : Global domain is distributed in smaller local domains.
[9019]	569	!!
[10425]	570	!! ** Action : - set for all knbi*knbj domains:
	571	!! kimppt : longitudinal index
	572	!! kjmppt : latitudinal index
	573	!! klci : first dimension
	574	!! klcj : second dimension
[9019]	575	!!----------------------------------------------------------------------
[13286]	576	INTEGER, INTENT(in ) :: kiglo, kjglo
	577	INTEGER, INTENT(in ) :: khls
[10425]	578	INTEGER, INTENT(in ) :: knbi, knbj
	579	INTEGER, INTENT( out) :: kimax, kjmax
	580	INTEGER, DIMENSION(knbi,knbj), OPTIONAL, INTENT( out) :: kimppt, kjmppt
	581	INTEGER, DIMENSION(knbi,knbj), OPTIONAL, INTENT( out) :: klci, klcj
	582	!
	583	INTEGER :: ji, jj
[14072]	584	INTEGER :: i2hls
[10542]	585	INTEGER :: iresti, irestj, irm, ijpjmin
[9019]	586	!!----------------------------------------------------------------------
[13286]	587	i2hls = 2*khls
[9019]	588	!
[10425]	589	#if defined key_nemocice_decomp
[13286]	590	kimax = ( nx_global+2-i2hls + (knbi-1) ) / knbi + i2hls ! first dim.
[14072]	591	kjmax = ( ny_global+2-i2hls + (knbj-1) ) / knbj + i2hls ! second dim.
[10425]	592	#else
[13286]	593	kimax = ( kiglo - i2hls + (knbi-1) ) / knbi + i2hls ! first dim.
	594	kjmax = ( kjglo - i2hls + (knbj-1) ) / knbj + i2hls ! second dim.
[10425]	595	#endif
	596	IF( .NOT. PRESENT(kimppt) ) RETURN
[9019]	597	!
[10425]	598	! 1. Dimension arrays for subdomains
	599	! -----------------------------------
	600	! Computation of local domain sizes klci() klcj()
[13286]	601	! These dimensions depend on global sizes knbi,knbj and kiglo,kjglo
[10425]	602	! The subdomains are squares lesser than or equal to the global
	603	! dimensions divided by the number of processors minus the overlap array.
[9019]	604	!
[13286]	605	iresti = 1 + MOD( kiglo - i2hls - 1 , knbi )
	606	irestj = 1 + MOD( kjglo - i2hls - 1 , knbj )
[9168]	607	!
[10425]	608	! Need to use kimax and kjmax here since jpi and jpj not yet defined
	609	#if defined key_nemocice_decomp
	610	! Change padding to be consistent with CICE
[13286]	611	klci(1:knbi-1,: ) = kimax
	612	klci( knbi ,: ) = kiglo - (knbi - 1) * (kimax - i2hls)
	613	klcj(: ,1:knbj-1) = kjmax
	614	klcj(: , knbj ) = kjglo - (knbj - 1) * (kjmax - i2hls)
[10425]	615	#else
	616	klci(1:iresti ,:) = kimax
	617	klci(iresti+1:knbi ,:) = kimax-1
[14848]	618	IF( MINVAL(klci) < 3*khls ) THEN
	619	WRITE(ctmp1,) ' mpp_basesplit: minimum value of jpi must be >= ', 3khls
[10542]	620	WRITE(ctmp2,*) ' We have ', MINVAL(klci)
[14848]	621	CALL ctl_stop( 'STOP', ctmp1, ctmp2 )
[10542]	622	ENDIF
[14433]	623	IF( l_NFold ) THEN
[10542]	624	! minimize the size of the last row to compensate for the north pole folding coast
[14433]	625	IF( c_NFtype == 'T' ) ijpjmin = 2+3*khls ! V and F folding must be outside of southern halos
	626	IF( c_NFtype == 'F' ) ijpjmin = 1+3*khls ! V and F folding must be outside of southern halos
	627	irm = knbj - irestj ! total number of lines to be removed
	628	klcj(:,knbj) = MAX( ijpjmin, kjmax-irm ) ! we must have jpj >= ijpjmin in the last row
	629	irm = irm - ( kjmax - klcj(1,knbj) ) ! remaining number of lines to remove
[13286]	630	irestj = knbj - 1 - irm
[10542]	631	klcj(:, irestj+1:knbj-1) = kjmax-1
	632	ELSE
[13286]	633	klcj(:, irestj+1:knbj ) = kjmax-1
[10542]	634	ENDIF
[13286]	635	klcj(:,1:irestj) = kjmax
[14848]	636	IF( MINVAL(klcj) < 3*khls ) THEN
	637	WRITE(ctmp1,) ' mpp_basesplit: minimum value of jpj must be >= ', 3khls
[10542]	638	WRITE(ctmp2,*) ' We have ', MINVAL(klcj)
	639	CALL ctl_stop( 'STOP', ctmp1, ctmp2 )
	640	ENDIF
[10425]	641	#endif
[3]	642
[10425]	643	! 2. Index arrays for subdomains
	644	! -------------------------------
	645	kimppt(:,:) = 1
	646	kjmppt(:,:) = 1
	647	!
	648	IF( knbi > 1 ) THEN
	649	DO jj = 1, knbj
	650	DO ji = 2, knbi
[13286]	651	kimppt(ji,jj) = kimppt(ji-1,jj) + klci(ji-1,jj) - i2hls
[10425]	652	END DO
	653	END DO
[9019]	654	ENDIF
	655	!
[10425]	656	IF( knbj > 1 )THEN
	657	DO jj = 2, knbj
	658	DO ji = 1, knbi
[13286]	659	kjmppt(ji,jj) = kjmppt(ji,jj-1) + klcj(ji,jj-1) - i2hls
[10425]	660	END DO
	661	END DO
	662	ENDIF
[14072]	663
[13286]	664	END SUBROUTINE mpp_basesplit
[9019]	665
	666
[13286]	667	SUBROUTINE bestpartition( knbij, knbi, knbj, knbcnt, ldlist )
[10425]	668	!!----------------------------------------------------------------------
[13286]	669	!! * ROUTINE bestpartition *
[10425]	670	!!
	671	!! ** Purpose :
	672	!!
	673	!! ** Method :
	674	!!----------------------------------------------------------------------
[13490]	675	INTEGER, INTENT(in ) :: knbij ! total number of subdomains (knbi*knbj)
[10425]	676	INTEGER, OPTIONAL, INTENT( out) :: knbi, knbj ! number if subdomains along i and j (knbi and knbj)
	677	INTEGER, OPTIONAL, INTENT( out) :: knbcnt ! number of land subdomains
	678	LOGICAL, OPTIONAL, INTENT(in ) :: ldlist ! .true.: print the list the best domain decompositions (with land)
	679	!
	680	INTEGER :: ji, jj, ii, iitarget
	681	INTEGER :: iszitst, iszjtst
	682	INTEGER :: isziref, iszjref
[13490]	683	INTEGER :: iszimin, iszjmin
[10425]	684	INTEGER :: inbij, iszij
[11536]	685	INTEGER :: inbimax, inbjmax, inbijmax, inbijold
[10425]	686	INTEGER :: isz0, isz1
	687	INTEGER, DIMENSION( :), ALLOCATABLE :: indexok
	688	INTEGER, DIMENSION( :), ALLOCATABLE :: inbi0, inbj0, inbij0 ! number of subdomains along i,j
	689	INTEGER, DIMENSION( :), ALLOCATABLE :: iszi0, iszj0, iszij0 ! max size of the subdomains along i,j
	690	INTEGER, DIMENSION( :), ALLOCATABLE :: inbi1, inbj1, inbij1 ! number of subdomains along i,j
	691	INTEGER, DIMENSION( :), ALLOCATABLE :: iszi1, iszj1, iszij1 ! max size of the subdomains along i,j
	692	LOGICAL :: llist
	693	LOGICAL, DIMENSION(:,:), ALLOCATABLE :: llmsk2d ! max size of the subdomains along i,j
[14433]	694	LOGICAL, DIMENSION(:,:), ALLOCATABLE :: llisOce ! - -
[10425]	695	REAL(wp):: zpropland
	696	!!----------------------------------------------------------------------
	697	!
	698	llist = .FALSE.
	699	IF( PRESENT(ldlist) ) llist = ldlist
	700
	701	CALL mpp_init_landprop( zpropland ) ! get the proportion of land point over the gloal domain
	702	inbij = NINT( REAL(knbij, wp) / ( 1.0 - zpropland ) ) ! define the largest possible value for jpni*jpnj
	703	!
	704	IF( llist ) THEN ; inbijmax = inbij*2
	705	ELSE ; inbijmax = inbij
	706	ENDIF
	707	!
	708	ALLOCATE(inbi0(inbijmax),inbj0(inbijmax),iszi0(inbijmax),iszj0(inbijmax))
	709	!
	710	inbimax = 0
	711	inbjmax = 0
[13490]	712	isziref = jpiglo*jpjglo+1 ! define a value that is larger than the largest possible
	713	iszjref = jpiglo*jpjglo+1
[10425]	714	!
[15296]	715	! WARNING, see also init_excl_landpt: minimum subdomain size defined here according to nn_hls (and not n_hlsmax)
	716	! --> If, one day, we want to use local halos largers than nn_hls, we must replace nn_hls by n_hlsmax
	717	!
[14848]	718	iszimin = 3*nn_hls ! minimum size of the MPI subdomain so halos are always adressing neighbor inner domain
	719	iszjmin = 3*nn_hls
[14433]	720	IF( c_NFtype == 'T' ) iszjmin = MAX(iszjmin, 2+3*nn_hls) ! V and F folding must be outside of southern halos
	721	IF( c_NFtype == 'F' ) iszjmin = MAX(iszjmin, 1+3*nn_hls) ! V and F folding must be outside of southern halos
[13490]	722	!
[10425]	723	! get the list of knbi that gives a smaller jpimax than knbi-1
	724	! get the list of knbj that gives a smaller jpjmax than knbj-1
[14072]	725	DO ji = 1, inbijmax
[10425]	726	#if defined key_nemocice_decomp
	727	iszitst = ( nx_global+2-2nn_hls + (ji-1) ) / ji + 2nn_hls ! first dim.
	728	#else
[13490]	729	iszitst = ( Ni0glo + (ji-1) ) / ji + 2*nn_hls ! max subdomain i-size
[10425]	730	#endif
[13490]	731	IF( iszitst < isziref .AND. iszitst >= iszimin ) THEN
[10425]	732	isziref = iszitst
	733	inbimax = inbimax + 1
	734	inbi0(inbimax) = ji
	735	iszi0(inbimax) = isziref
	736	ENDIF
	737	#if defined key_nemocice_decomp
	738	iszjtst = ( ny_global+2-2nn_hls + (ji-1) ) / ji + 2nn_hls ! first dim.
	739	#else
[13490]	740	iszjtst = ( Nj0glo + (ji-1) ) / ji + 2*nn_hls ! max subdomain j-size
[10425]	741	#endif
[13490]	742	IF( iszjtst < iszjref .AND. iszjtst >= iszjmin ) THEN
[10425]	743	iszjref = iszjtst
	744	inbjmax = inbjmax + 1
	745	inbj0(inbjmax) = ji
	746	iszj0(inbjmax) = iszjref
	747	ENDIF
	748	END DO
[14848]	749	IF( inbimax == 0 ) THEN
	750	WRITE(ctmp1,'(a,i2,a,i2)') ' mpp_ini bestpartition: Ni0glo (', Ni0glo, ') is too small to be used with nn_hls = ', nn_hls
	751	CALL ctl_stop( 'STOP', ctmp1 )
	752	ENDIF
	753	IF( inbjmax == 0 ) THEN
	754	WRITE(ctmp1,'(a,i2,a,i2)') ' mpp_ini bestpartition: Nj0glo (', Nj0glo, ') is too small to be used with nn_hls = ', nn_hls
	755	CALL ctl_stop( 'STOP', ctmp1 )
	756	ENDIF
[10425]	757
	758	! combine these 2 lists to get all possible knbi*knbj < inbijmax
	759	ALLOCATE( llmsk2d(inbimax,inbjmax) )
	760	DO jj = 1, inbjmax
	761	DO ji = 1, inbimax
	762	IF ( inbi0(ji) * inbj0(jj) <= inbijmax ) THEN ; llmsk2d(ji,jj) = .TRUE.
	763	ELSE ; llmsk2d(ji,jj) = .FALSE.
	764	ENDIF
	765	END DO
	766	END DO
	767	isz1 = COUNT(llmsk2d)
	768	ALLOCATE( inbi1(isz1), inbj1(isz1), iszi1(isz1), iszj1(isz1) )
	769	ii = 0
	770	DO jj = 1, inbjmax
	771	DO ji = 1, inbimax
	772	IF( llmsk2d(ji,jj) .EQV. .TRUE. ) THEN
	773	ii = ii + 1
	774	inbi1(ii) = inbi0(ji)
	775	inbj1(ii) = inbj0(jj)
	776	iszi1(ii) = iszi0(ji)
	777	iszj1(ii) = iszj0(jj)
[14433]	778	ENDIF
[10425]	779	END DO
	780	END DO
	781	DEALLOCATE( inbi0, inbj0, iszi0, iszj0 )
	782	DEALLOCATE( llmsk2d )
	783
	784	ALLOCATE( inbij1(isz1), iszij1(isz1) )
	785	inbij1(:) = inbi1(:) * inbj1(:)
	786	iszij1(:) = iszi1(:) * iszj1(:)
	787
[13286]	788	! if there is no land and no print
[10425]	789	IF( .NOT. llist .AND. numbot == -1 .AND. numbdy == -1 ) THEN
	790	! get the smaller partition which gives the smallest subdomain size
	791	ii = MINLOC(inbij1, mask = iszij1 == MINVAL(iszij1), dim = 1)
	792	knbi = inbi1(ii)
	793	knbj = inbj1(ii)
	794	IF(PRESENT(knbcnt)) knbcnt = 0
	795	DEALLOCATE( inbi1, inbj1, inbij1, iszi1, iszj1, iszij1 )
	796	RETURN
	797	ENDIF
	798
	799	! extract only the partitions which reduce the subdomain size in comparison with smaller partitions
	800	ALLOCATE( indexok(isz1) ) ! to store indices of the best partitions
[14072]	801	isz0 = 0 ! number of best partitions
[10425]	802	inbij = 1 ! start with the min value of inbij1 => 1
[13490]	803	iszij = jpiglo*jpjglo+1 ! default: larger than global domain
[10425]	804	DO WHILE( inbij <= inbijmax ) ! if we did not reach the max of inbij1
	805	ii = MINLOC(iszij1, mask = inbij1 == inbij, dim = 1) ! warning: send back the first occurence if multiple results
	806	IF ( iszij1(ii) < iszij ) THEN
[13490]	807	ii = MINLOC( iszi1+iszj1, mask = iszij1 == iszij1(ii) .AND. inbij1 == inbij, dim = 1) ! select the smaller perimeter if multiple min
[10425]	808	isz0 = isz0 + 1
	809	indexok(isz0) = ii
	810	iszij = iszij1(ii)
	811	ENDIF
	812	inbij = MINVAL(inbij1, mask = inbij1 > inbij) ! warning: return largest integer value if mask = .false. everywhere
	813	END DO
	814	DEALLOCATE( inbij1, iszij1 )
	815
	816	! keep only the best partitions (sorted by increasing order of subdomains number and decreassing subdomain size)
	817	ALLOCATE( inbi0(isz0), inbj0(isz0), iszi0(isz0), iszj0(isz0) )
	818	DO ji = 1, isz0
	819	ii = indexok(ji)
	820	inbi0(ji) = inbi1(ii)
	821	inbj0(ji) = inbj1(ii)
	822	iszi0(ji) = iszi1(ii)
	823	iszj0(ji) = iszj1(ii)
	824	END DO
	825	DEALLOCATE( indexok, inbi1, inbj1, iszi1, iszj1 )
	826
[11536]	827	IF( llist ) THEN
[10425]	828	IF(lwp) THEN
	829	WRITE(numout,*)
[11536]	830	WRITE(numout,*) ' For your information:'
	831	WRITE(numout,*) ' list of the best partitions including land supression'
	832	WRITE(numout,*) ' -----------------------------------------------------'
[10425]	833	WRITE(numout,*)
[14433]	834	ENDIF
[11536]	835	ji = isz0 ! initialization with the largest value
[14433]	836	ALLOCATE( llisOce(inbi0(ji), inbj0(ji)) )
	837	CALL mpp_is_ocean( llisOce ) ! Warning: must be call by all cores (call mpp_sum)
	838	inbijold = COUNT(llisOce)
	839	DEALLOCATE( llisOce )
[11536]	840	DO ji =isz0-1,1,-1
[14433]	841	ALLOCATE( llisOce(inbi0(ji), inbj0(ji)) )
	842	CALL mpp_is_ocean( llisOce ) ! Warning: must be call by all cores (call mpp_sum)
	843	inbij = COUNT(llisOce)
	844	DEALLOCATE( llisOce )
[11536]	845	IF(lwp .AND. inbij < inbijold) THEN
	846	WRITE(numout,'(a, i6, a, i6, a, f4.1, a, i9, a, i6, a, i6, a)') &
	847	& 'nb_cores oce: ', inbij, ', land domains excluded: ', inbi0(ji)*inbj0(ji) - inbij, &
	848	& ' (', REAL(inbi0(ji)inbj0(ji) - inbij,wp) / REAL(inbi0(ji)inbj0(ji),wp) *100., &
	849	& '%), largest oce domain: ', iszi0(ji)*iszj0(ji), ' ( ', iszi0(ji),' x ', iszj0(ji), ' )'
	850	inbijold = inbij
[14433]	851	ENDIF
[10425]	852	END DO
	853	DEALLOCATE( inbi0, inbj0, iszi0, iszj0 )
[11536]	854	IF(lwp) THEN
	855	WRITE(numout,*)
	856	WRITE(numout,*) ' -----------------------------------------------------------'
	857	ENDIF
	858	CALL mppsync
	859	CALL mppstop( ld_abort = .TRUE. )
[10425]	860	ENDIF
[14072]	861
[10425]	862	DEALLOCATE( iszi0, iszj0 )
	863	inbij = inbijmax + 1 ! default: larger than possible
	864	ii = isz0+1 ! start from the end of the list (smaller subdomains)
	865	DO WHILE( inbij > knbij ) ! while the number of ocean subdomains exceed the number of procs
[14072]	866	ii = ii -1
[14433]	867	ALLOCATE( llisOce(inbi0(ii), inbj0(ii)) )
	868	CALL mpp_is_ocean( llisOce ) ! must be done by all core
	869	inbij = COUNT(llisOce)
	870	DEALLOCATE( llisOce )
[10425]	871	END DO
	872	knbi = inbi0(ii)
	873	knbj = inbj0(ii)
	874	IF(PRESENT(knbcnt)) knbcnt = knbi * knbj - inbij
	875	DEALLOCATE( inbi0, inbj0 )
	876	!
[13286]	877	END SUBROUTINE bestpartition
[14072]	878
	879
[10425]	880	SUBROUTINE mpp_init_landprop( propland )
	881	!!----------------------------------------------------------------------
	882	!! * ROUTINE mpp_init_landprop *
	883	!!
	884	!! ** Purpose : the the proportion of land points in the surface land-sea mask
	885	!!
[13286]	886	!! ** Method : read iproc strips (of length Ni0glo) of the land-sea mask
[10425]	887	!!----------------------------------------------------------------------
[10436]	888	REAL(wp), INTENT( out) :: propland ! proportion of land points in the global domain (between 0 and 1)
[10425]	889	!
	890	INTEGER, DIMENSION(jpni*jpnj) :: kusedom_1d
[10436]	891	INTEGER :: inboce, iarea
[10425]	892	INTEGER :: iproc, idiv, ijsz
[10436]	893	INTEGER :: ijstr
[10425]	894	LOGICAL, ALLOCATABLE, DIMENSION(:,:) :: lloce
	895	!!----------------------------------------------------------------------
	896	! do nothing if there is no land-sea mask
	897	IF( numbot == -1 .and. numbdy == -1 ) THEN
	898	propland = 0.
	899	RETURN
	900	ENDIF
	901
[14072]	902	! number of processes reading the bathymetry file
[13286]	903	iproc = MINVAL( (/mppsize, Nj0glo/2, 100/) ) ! read a least 2 lines, no more that 100 processes reading at the same time
[14072]	904
[10436]	905	! we want to read iproc strips of the land-sea mask. -> pick up iproc processes every idiv processes starting at 1
[10425]	906	IF( iproc == 1 ) THEN ; idiv = mppsize
	907	ELSE ; idiv = ( mppsize - 1 ) / ( iproc - 1 )
	908	ENDIF
[10436]	909
	910	iarea = (narea-1)/idiv ! involed process number (starting counting at 0)
	911	IF( MOD( narea-1, idiv ) == 0 .AND. iarea < iproc ) THEN ! beware idiv can be = to 1
[10425]	912	!
[13286]	913	ijsz = Nj0glo / iproc ! width of the stripe to read
	914	IF( iarea < MOD(Nj0glo,iproc) ) ijsz = ijsz + 1
	915	ijstr = iarea*(Nj0glo/iproc) + MIN(iarea, MOD(Nj0glo,iproc)) + 1 ! starting j position of the reading
[10425]	916	!
[13286]	917	ALLOCATE( lloce(Ni0glo, ijsz) ) ! allocate the strip
[14433]	918	CALL read_mask( 1, ijstr, Ni0glo, ijsz, lloce )
[10436]	919	inboce = COUNT(lloce) ! number of ocean point in the stripe
[10425]	920	DEALLOCATE(lloce)
	921	!
	922	ELSE
	923	inboce = 0
	924	ENDIF
[10436]	925	CALL mpp_sum( 'mppini', inboce ) ! total number of ocean points over the global domain
[10425]	926	!
[14072]	927	propland = REAL( Ni0gloNj0glo - inboce, wp ) / REAL( Ni0gloNj0glo, wp )
[10425]	928	!
	929	END SUBROUTINE mpp_init_landprop
[14072]	930
	931
[14433]	932	SUBROUTINE mpp_is_ocean( ldIsOce )
[10425]	933	!!----------------------------------------------------------------------
[13286]	934	!! * ROUTINE mpp_is_ocean *
[10425]	935	!!
[13286]	936	!! ** Purpose : Check for a mpi domain decomposition inbi x inbj which
	937	!! subdomains, including 1 halo (even if nn_hls>1), contain
	938	!! at least 1 ocean point.
	939	!! We must indeed ensure that each subdomain that is a neighbour
[14433]	940	!! of a land subdomain, has only land points on its boundary
[13286]	941	!! (inside the inner subdomain) with the land subdomain.
	942	!! This is needed to get the proper bondary conditions on
	943	!! a subdomain with a closed boundary.
[10425]	944	!!
[13286]	945	!! ** Method : read inbj strips (of length Ni0glo) of the land-sea mask
[10425]	946	!!----------------------------------------------------------------------
[14433]	947	LOGICAL, DIMENSION(:,:), INTENT( out) :: ldIsOce ! .true. if a sub domain constains 1 ocean point
[10425]	948	!
[10436]	949	INTEGER :: idiv, iimax, ijmax, iarea
[13286]	950	INTEGER :: inbi, inbj, inx, iny, inry, isty
[10560]	951	INTEGER :: ji, jn
[13286]	952	INTEGER, ALLOCATABLE, DIMENSION(:,:) :: inboce ! number oce oce pint in each mpi subdomain
	953	INTEGER, ALLOCATABLE, DIMENSION(: ) :: inboce_1d
	954	INTEGER, ALLOCATABLE, DIMENSION(:,:) :: iimppt, ijpi
	955	INTEGER, ALLOCATABLE, DIMENSION(:,:) :: ijmppt, ijpj
[14072]	956	LOGICAL, ALLOCATABLE, DIMENSION(:,:) :: lloce ! lloce(i,j) = .true. if the point (i,j) is ocean
[10425]	957	!!----------------------------------------------------------------------
	958	! do nothing if there is no land-sea mask
	959	IF( numbot == -1 .AND. numbdy == -1 ) THEN
[14433]	960	ldIsOce(:,:) = .TRUE.
[10425]	961	RETURN
	962	ENDIF
[13286]	963	!
[14433]	964	inbi = SIZE( ldIsOce, dim = 1 )
	965	inbj = SIZE( ldIsOce, dim = 2 )
[13286]	966	!
	967	! we want to read inbj strips of the land-sea mask. -> pick up inbj processes every idiv processes starting at 1
	968	IF ( inbj == 1 ) THEN ; idiv = mppsize
	969	ELSE IF ( mppsize < inbj ) THEN ; idiv = 1
	970	ELSE ; idiv = ( mppsize - 1 ) / ( inbj - 1 )
[10425]	971	ENDIF
[13286]	972	!
	973	ALLOCATE( inboce(inbi,inbj), inboce_1d(inbi*inbj) )
[10436]	974	inboce(:,:) = 0 ! default no ocean point found
[13286]	975	!
	976	DO jn = 0, (inbj-1)/mppsize ! if mppsize < inbj : more strips than mpi processes (because of potential land domains)
[10425]	977	!
[13286]	978	iarea = (narea-1)/idiv + jn * mppsize + 1 ! involed process number (starting counting at 1)
	979	IF( MOD( narea-1, idiv ) == 0 .AND. iarea <= inbj ) THEN ! beware idiv can be = to 1
[10560]	980	!
[13286]	981	ALLOCATE( iimppt(inbi,inbj), ijmppt(inbi,inbj), ijpi(inbi,inbj), ijpj(inbi,inbj) )
	982	CALL mpp_basesplit( Ni0glo, Nj0glo, 0, inbi, inbj, iimax, ijmax, iimppt, ijmppt, ijpi, ijpj )
[10560]	983	!
[13286]	984	inx = Ni0glo + 2 ; iny = ijpj(1,iarea) + 2 ! strip size + 1 halo on each direction (even if nn_hls>1)
	985	ALLOCATE( lloce(inx, iny) ) ! allocate the strip
	986	inry = iny - COUNT( (/ iarea == 1, iarea == inbj /) ) ! number of point to read in y-direction
	987	isty = 1 + COUNT( (/ iarea == 1 /) ) ! read from the first or the second line?
[14433]	988	CALL read_mask( 1, ijmppt(1,iarea) - 2 + isty, Ni0glo, inry, lloce(2:inx-1, isty:inry+isty-1) ) ! read the strip
[14072]	989	!
[13286]	990	IF( iarea == 1 ) THEN ! the first line was not read
[14433]	991	IF( l_Jperio ) THEN ! north-south periodocity
	992	CALL read_mask( 1, Nj0glo, Ni0glo, 1, lloce(2:inx-1, 1) ) ! read the last line -> first line of lloce
[13286]	993	ELSE
	994	lloce(2:inx-1, 1) = .FALSE. ! closed boundary
	995	ENDIF
	996	ENDIF
	997	IF( iarea == inbj ) THEN ! the last line was not read
[14433]	998	IF( l_Jperio ) THEN ! north-south periodocity
	999	CALL read_mask( 1, 1, Ni0glo, 1, lloce(2:inx-1,iny) ) ! read the first line -> last line of lloce
	1000	ELSEIF( c_NFtype == 'T' ) THEN ! north-pole folding T-pivot, T-point
[13286]	1001	lloce(2,iny) = lloce(2,iny-2) ! here we have 1 halo (even if nn_hls>1)
	1002	DO ji = 3,inx-1
	1003	lloce(ji,iny ) = lloce(inx-ji+2,iny-2) ! ok, we have at least 3 lines
	1004	END DO
	1005	DO ji = inx/2+2,inx-1
	1006	lloce(ji,iny-1) = lloce(inx-ji+2,iny-1)
	1007	END DO
[14433]	1008	ELSEIF( c_NFtype == 'F' ) THEN ! north-pole folding F-pivot, T-point, 1 halo
[13286]	1009	lloce(inx/2+1,iny-1) = lloce(inx/2,iny-1) ! here we have 1 halo (even if nn_hls>1)
	1010	lloce(inx -1,iny-1) = lloce(2 ,iny-1)
	1011	DO ji = 2,inx-1
	1012	lloce(ji,iny) = lloce(inx-ji+1,iny-1)
	1013	END DO
	1014	ELSE ! closed boundary
	1015	lloce(2:inx-1,iny) = .FALSE.
	1016	ENDIF
	1017	ENDIF
	1018	! ! first and last column were not read
[14433]	1019	IF( l_Iperio ) THEN
[13286]	1020	lloce(1,:) = lloce(inx-1,:) ; lloce(inx,:) = lloce(2,:) ! east-west periodocity
	1021	ELSE
	1022	lloce(1,:) = .FALSE. ; lloce(inx,:) = .FALSE. ! closed boundary
	1023	ENDIF
	1024	!
	1025	DO ji = 1, inbi
	1026	inboce(ji,iarea) = COUNT( lloce(iimppt(ji,1):iimppt(ji,1)+ijpi(ji,1)+1,:) ) ! lloce as 2 points more than Ni0glo
[10560]	1027	END DO
	1028	!
	1029	DEALLOCATE(lloce)
[13286]	1030	DEALLOCATE(iimppt, ijmppt, ijpi, ijpj)
[10560]	1031	!
	1032	ENDIF
	1033	END DO
[14072]	1034
[13286]	1035	inboce_1d = RESHAPE(inboce, (/ inbi*inbj /))
[10425]	1036	CALL mpp_sum( 'mppini', inboce_1d )
[13286]	1037	inboce = RESHAPE(inboce_1d, (/inbi, inbj/))
[14433]	1038	ldIsOce(:,:) = inboce(:,:) /= 0
[13286]	1039	DEALLOCATE(inboce, inboce_1d)
[10425]	1040	!
[13286]	1041	END SUBROUTINE mpp_is_ocean
[14072]	1042
	1043
[14433]	1044	SUBROUTINE read_mask( kistr, kjstr, kicnt, kjcnt, ldoce )
[10425]	1045	!!----------------------------------------------------------------------
[14433]	1046	!! * ROUTINE read_mask *
[10425]	1047	!!
	1048	!! ** Purpose : Read relevant bathymetric information in order to
	1049	!! provide a land/sea mask used for the elimination
	1050	!! of land domains, in an mpp computation.
	1051	!!
[13286]	1052	!! ** Method : read stipe of size (Ni0glo,...)
[10425]	1053	!!----------------------------------------------------------------------
[14433]	1054	INTEGER , INTENT(in ) :: kistr, kjstr ! starting i and j position of the reading
	1055	INTEGER , INTENT(in ) :: kicnt, kjcnt ! number of points to read in i and j directions
	1056	LOGICAL, DIMENSION(kicnt,kjcnt), INTENT( out) :: ldoce ! ldoce(i,j) = .true. if the point (i,j) is ocean
[10425]	1057	!
[14433]	1058	INTEGER :: inumsave ! local logical unit
	1059	REAL(wp), DIMENSION(kicnt,kjcnt) :: zbot, zbdy
[10425]	1060	!!----------------------------------------------------------------------
	1061	!
	1062	inumsave = numout ; numout = numnul ! redirect all print to /dev/null
	1063	!
[14072]	1064	IF( numbot /= -1 ) THEN
[14433]	1065	CALL iom_get( numbot, jpdom_unknown, 'bottom_level', zbot, kstart = (/kistr,kjstr/), kcount = (/kicnt, kjcnt/) )
[10425]	1066	ELSE
[13286]	1067	zbot(:,:) = 1._wp ! put a non-null value
[10425]	1068	ENDIF
[13286]	1069	!
[14072]	1070	IF( numbdy /= -1 ) THEN ! Adjust with bdy_msk if it exists
[14433]	1071	CALL iom_get ( numbdy, jpdom_unknown, 'bdy_msk', zbdy, kstart = (/kistr,kjstr/), kcount = (/kicnt, kjcnt/) )
[10425]	1072	zbot(:,:) = zbot(:,:) * zbdy(:,:)
	1073	ENDIF
	1074	!
[14433]	1075	ldoce(:,:) = NINT(zbot(:,:)) > 0
[10425]	1076	numout = inumsave
	1077	!
[14433]	1078	END SUBROUTINE read_mask
[10425]	1079
	1080
[14433]	1081	SUBROUTINE mpp_getnum( ldIsOce, kproc, kipos, kjpos )
[88]	1082	!!----------------------------------------------------------------------
[13286]	1083	!! * ROUTINE mpp_getnum *
[88]	1084	!!
[13286]	1085	!! ** Purpose : give a number to each MPI subdomains (starting at 0)
	1086	!!
	1087	!! ** Method : start from bottom left. First skip land subdomain, and finally use them if needed
	1088	!!----------------------------------------------------------------------
[14433]	1089	LOGICAL, DIMENSION(:,:), INTENT(in ) :: ldIsOce ! F if land process
	1090	INTEGER, DIMENSION(:,:), INTENT( out) :: kproc ! subdomain number (-1 if not existing, starting at 0)
[13286]	1091	INTEGER, DIMENSION( :), INTENT( out) :: kipos ! i-position of the subdomain (from 1 to jpni)
	1092	INTEGER, DIMENSION( :), INTENT( out) :: kjpos ! j-position of the subdomain (from 1 to jpnj)
	1093	!
	1094	INTEGER :: ii, ij, jarea, iarea0
	1095	INTEGER :: icont, i2add , ini, inj, inij
	1096	!!----------------------------------------------------------------------
	1097	!
[14433]	1098	ini = SIZE(ldIsOce, dim = 1)
	1099	inj = SIZE(ldIsOce, dim = 2)
[13286]	1100	inij = SIZE(kipos)
	1101	!
	1102	! specify which subdomains are oce subdomains; other are land subdomains
	1103	kproc(:,:) = -1
	1104	icont = -1
	1105	DO jarea = 1, ini*inj
	1106	iarea0 = jarea - 1
	1107	ii = 1 + MOD(iarea0,ini)
	1108	ij = 1 + iarea0/ini
[14433]	1109	IF( ldIsOce(ii,ij) ) THEN
[13286]	1110	icont = icont + 1
	1111	kproc(ii,ij) = icont
	1112	kipos(icont+1) = ii
	1113	kjpos(icont+1) = ij
	1114	ENDIF
	1115	END DO
	1116	! if needed add some land subdomains to reach inij active subdomains
[14433]	1117	i2add = inij - COUNT( ldIsOce )
[13286]	1118	DO jarea = 1, ini*inj
	1119	iarea0 = jarea - 1
	1120	ii = 1 + MOD(iarea0,ini)
	1121	ij = 1 + iarea0/ini
[14433]	1122	IF( .NOT. ldIsOce(ii,ij) .AND. i2add > 0 ) THEN
[13286]	1123	icont = icont + 1
	1124	kproc(ii,ij) = icont
	1125	kipos(icont+1) = ii
	1126	kjpos(icont+1) = ij
	1127	i2add = i2add - 1
	1128	ENDIF
	1129	END DO
	1130	!
	1131	END SUBROUTINE mpp_getnum
	1132
	1133
[14433]	1134	SUBROUTINE init_excl_landpt
	1135	!!----------------------------------------------------------------------
	1136	!! * ROUTINE *
	1137	!!
	1138	!! ** Purpose : exclude exchanges which contain only land points
	1139	!!
	1140	!! ** Method : if a send or receive buffer constains only land point we
	1141	!! flag off the corresponding communication
	1142	!! Warning: this selection depend on the halo size -> loop on halo size
	1143	!!
	1144	!!----------------------------------------------------------------------
	1145	INTEGER :: inumsave
	1146	INTEGER :: jh
	1147	INTEGER :: ipi, ipj
	1148	INTEGER :: iiwe, iiea, iist, iisz
	1149	INTEGER :: ijso, ijno, ijst, ijsz
	1150	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zmsk
	1151	LOGICAL , DIMENSION(Ni_0,Nj_0,1) :: lloce
	1152	!!----------------------------------------------------------------------
	1153	!
	1154	! read the land-sea mask on the inner domain
	1155	CALL read_mask( nimpp, njmpp, Ni_0, Nj_0, lloce(:,:,1) )
	1156	!
	1157	! Here we look only at communications excluding the NP folding.
[15267]	1158	! --> we switch off lbcnfd at this stage (init_nfdcom called after init_excl_landpt)...
[14433]	1159	l_IdoNFold = .FALSE.
	1160	!
[15296]	1161	! WARNING, see also bestpartition: minimum subdomain size defined in bestpartition according to nn_hls.
	1162	! If, one day, we want to use local halos largers than nn_hls, we must replace nn_hls by n_hlsmax in bestpartition
	1163	!
	1164	DO jh = 1, MIN(nn_hls, n_hlsmax) ! different halo size
[14433]	1165	!
	1166	ipi = Ni_0 + 2*jh ! local domain size
	1167	ipj = Nj_0 + 2*jh
	1168	!
	1169	ALLOCATE( zmsk(ipi,ipj) )
	1170	zmsk(jh+1:jh+Ni_0,jh+1:jh+Nj_0) = REAL(COUNT(lloce, dim = 3), wp) ! define inner domain -> need REAL to use lbclnk
	1171	CALL lbc_lnk('mppini', zmsk, 'T', 1._wp, khls = jh) ! fill halos
[15299]	1172	! Beware, coastal F points can be used in the code -> we may need communications for these points F points even if tmask = 0
	1173	! -> the mask we must use here is equal to 1 as soon as one of the 4 neighbours is oce (sum of the mask, not multiplication)
	1174	zmsk(jh+1:jh+Ni_0,jh+1:jh+Nj_0) = zmsk(jh+1:jh+Ni_0,jh+1 :jh+Nj_0 ) + zmsk(jh+1+1:jh+Ni_0+1,jh+1 :jh+Nj_0 ) &
	1175	& + zmsk(jh+1:jh+Ni_0,jh+1+1:jh+Nj_0+1) + zmsk(jh+1+1:jh+Ni_0+1,jh+1+1:jh+Nj_0+1)
	1176	CALL lbc_lnk('mppini', zmsk, 'T', 1._wp, khls = jh) ! fill halos again!
[14433]	1177	!
[15296]	1178	iiwe = jh ; iiea = Ni_0 ! bottom-left corner - 1 of the sent data
[14433]	1179	ijso = jh ; ijno = Nj_0
	1180	IF( nn_comm == 1 ) THEN
	1181	iist = 0 ; iisz = ipi
[15302]	1182	ijst = jh ; ijsz = Nj_0
[14433]	1183	ELSE
	1184	iist = jh ; iisz = Ni_0
	1185	ijst = jh ; ijsz = Nj_0
	1186	ENDIF
	1187	IF( nn_comm == 1 ) THEN ! SM: NOT WORKING FOR NEIGHBOURHOOD COLLECTIVE COMMUNICATIONS, I DON'T KNOW WHY...
	1188	! do not send if we send only land points
	1189	IF( NINT(SUM( zmsk(iiwe+1:iiwe+jh ,ijst+1:ijst+ijsz) )) == 0 ) mpiSnei(jh,jpwe) = -1
	1190	IF( NINT(SUM( zmsk(iiea+1:iiea+jh ,ijst+1:ijst+ijsz) )) == 0 ) mpiSnei(jh,jpea) = -1
	1191	IF( NINT(SUM( zmsk(iist+1:iist+iisz,ijso+1:ijso+jh ) )) == 0 ) mpiSnei(jh,jpso) = -1
	1192	IF( NINT(SUM( zmsk(iist+1:iist+iisz,ijno+1:ijno+jh ) )) == 0 ) mpiSnei(jh,jpno) = -1
	1193	IF( NINT(SUM( zmsk(iiwe+1:iiwe+jh ,ijso+1:ijso+jh ) )) == 0 ) mpiSnei(jh,jpsw) = -1
	1194	IF( NINT(SUM( zmsk(iiea+1:iiea+jh ,ijso+1:ijso+jh ) )) == 0 ) mpiSnei(jh,jpse) = -1
	1195	IF( NINT(SUM( zmsk(iiwe+1:iiwe+jh ,ijno+1:ijno+jh ) )) == 0 ) mpiSnei(jh,jpnw) = -1
	1196	IF( NINT(SUM( zmsk(iiea+1:iiea+jh ,ijno+1:ijno+jh ) )) == 0 ) mpiSnei(jh,jpne) = -1
	1197	!
[15296]	1198	iiwe = iiwe-jh ; iiea = iiea+jh ! bottom-left corner - 1 of the received data
[14433]	1199	ijso = ijso-jh ; ijno = ijno+jh
	1200	! do not send if we send only land points
	1201	IF( NINT(SUM( zmsk(iiwe+1:iiwe+jh ,ijst+1:ijst+ijsz) )) == 0 ) mpiRnei(jh,jpwe) = -1
	1202	IF( NINT(SUM( zmsk(iiea+1:iiea+jh ,ijst+1:ijst+ijsz) )) == 0 ) mpiRnei(jh,jpea) = -1
	1203	IF( NINT(SUM( zmsk(iist+1:iist+iisz,ijso+1:ijso+jh ) )) == 0 ) mpiRnei(jh,jpso) = -1
	1204	IF( NINT(SUM( zmsk(iist+1:iist+iisz,ijno+1:ijno+jh ) )) == 0 ) mpiRnei(jh,jpno) = -1
	1205	IF( NINT(SUM( zmsk(iiwe+1:iiwe+jh ,ijso+1:ijso+jh ) )) == 0 ) mpiRnei(jh,jpsw) = -1
	1206	IF( NINT(SUM( zmsk(iiea+1:iiea+jh ,ijso+1:ijso+jh ) )) == 0 ) mpiRnei(jh,jpse) = -1
	1207	IF( NINT(SUM( zmsk(iiwe+1:iiwe+jh ,ijno+1:ijno+jh ) )) == 0 ) mpiRnei(jh,jpnw) = -1
	1208	IF( NINT(SUM( zmsk(iiea+1:iiea+jh ,ijno+1:ijno+jh ) )) == 0 ) mpiRnei(jh,jpne) = -1
	1209	ENDIF
	1210	!
	1211	! Specific (and rare) problem in corner treatment because we do 1st West-East comm, next South-North comm
	1212	IF( nn_comm == 1 ) THEN
	1213	IF( mpiSnei(jh,jpwe) > -1 ) mpiSnei(jh, (/jpsw,jpnw/) ) = -1 ! SW and NW corners already sent through West nei
	1214	IF( mpiSnei(jh,jpea) > -1 ) mpiSnei(jh, (/jpse,jpne/) ) = -1 ! SE and NE corners already sent through East nei
	1215	IF( mpiRnei(jh,jpso) > -1 ) mpiRnei(jh, (/jpsw,jpse/) ) = -1 ! SW and SE corners will be received through South nei
	1216	IF( mpiRnei(jh,jpno) > -1 ) mpiRnei(jh, (/jpnw,jpne/) ) = -1 ! NW and NE corners will be received through North nei
	1217	ENDIF
	1218	!
	1219	DEALLOCATE( zmsk )
	1220	!
	1221	CALL mpp_ini_nc(jh) ! Initialize/Update communicator for neighbourhood collective communications
	1222	!
	1223	END DO
	1224
	1225	END SUBROUTINE init_excl_landpt
	1226
	1227
[13286]	1228	SUBROUTINE init_ioipsl
	1229	!!----------------------------------------------------------------------
	1230	!! * ROUTINE init_ioipsl *
	1231	!!
[14072]	1232	!! ** Purpose :
[88]	1233	!!
[14072]	1234	!! ** Method :
[88]	1235	!!
	1236	!! History :
[14072]	1237	!! 9.0 ! 04-03 (G. Madec ) MPP-IOIPSL
[1238]	1238	!! " " ! 08-12 (A. Coward) addition in case of jpni*jpnj < jpnij
[88]	1239	!!----------------------------------------------------------------------
[2715]	1240	INTEGER, DIMENSION(2) :: iglo, iloc, iabsf, iabsl, ihals, ihale, idid
[88]	1241	!!----------------------------------------------------------------------
[352]	1242
[1238]	1243	! The domain is split only horizontally along i- or/and j- direction
	1244	! So we need at the most only 1D arrays with 2 elements.
	1245	! Set idompar values equivalent to the jpdom_local_noextra definition
	1246	! used in IOM. This works even if jpnij .ne. jpni*jpnj.
[13286]	1247	iglo( :) = (/ Ni0glo, Nj0glo /)
	1248	iloc( :) = (/ Ni_0 , Nj_0 /)
	1249	iabsf(:) = (/ Nis0 , Njs0 /) + (/ nimpp, njmpp /) - 1 - nn_hls ! corresponds to mig0(Nis0) but mig0 is not yet defined!
[88]	1250	iabsl(:) = iabsf(:) + iloc(:) - 1
[13286]	1251	ihals(:) = (/ 0 , 0 /)
	1252	ihale(:) = (/ 0 , 0 /)
	1253	idid( :) = (/ 1 , 2 /)
[352]	1254
[88]	1255	IF(lwp) THEN
[516]	1256	WRITE(numout,*)
[13286]	1257	WRITE(numout,*) 'mpp init_ioipsl : iloc = ', iloc
	1258	WRITE(numout,*) '~~~~~~~~~~~~~~~ iabsf = ', iabsf
	1259	WRITE(numout,*) ' ihals = ', ihals
	1260	WRITE(numout,*) ' ihale = ', ihale
[88]	1261	ENDIF
[2715]	1262	!
[14275]	1263	CALL flio_dom_set ( jpnij, narea-1, idid, iglo, iloc, iabsf, iabsl, ihals, ihale, 'BOX', nidom)
[2715]	1264	!
[14072]	1265	END SUBROUTINE init_ioipsl
[88]	1266
[9436]	1267
[15267]	1268	SUBROUTINE init_nfdcom( ldwrtlay, knum )
[9436]	1269	!!----------------------------------------------------------------------
[13286]	1270	!! * ROUTINE init_nfdcom *
[14072]	1271	!! ** Purpose : Setup for north fold exchanges with explicit
[9436]	1272	!! point-to-point messaging
	1273	!!
	1274	!! ** Method : Initialization of the northern neighbours lists.
	1275	!!----------------------------------------------------------------------
	1276	!! 1.0 ! 2011-10 (A. C. Coward, NOCS & J. Donners, PRACE)
[14072]	1277	!! 2.0 ! 2013-06 Setup avoiding MPI communication (I. Epicoco, S. Mocavero, CMCC)
[15267]	1278	!! 3.0 ! 2021-09 complete rewrite using informations from gather north fold
[9436]	1279	!!----------------------------------------------------------------------
[15267]	1280	LOGICAL, INTENT(in ) :: ldwrtlay ! true if additional prints in layout.dat
	1281	INTEGER, INTENT(in ) :: knum ! layout.dat unit
	1282	!
	1283	REAL(wp), DIMENSION(jpi,jpj,2,4) :: zinfo
	1284	INTEGER , DIMENSION(10) :: irknei ! too many elements but safe...
	1285	INTEGER :: ji, jj, jg, jn ! dummy loop indices
	1286	LOGICAL :: lnew
[9436]	1287	!!----------------------------------------------------------------------
	1288	!
[15267]	1289	IF (lwp) THEN
	1290	WRITE(numout,*)
	1291	WRITE(numout,*) ' ==>>> North fold boundary prepared for jpni >1'
	1292	ENDIF
[9436]	1293	!
[15267]	1294	CALL mpp_ini_northgather ! we need to init the nfd with gathering in all cases as it is used to define the no-gather case
[14433]	1295	!
[15267]	1296	IF(ldwrtlay) THEN ! additional prints in layout.dat
	1297	WRITE(knum,*)
	1298	WRITE(knum,*)
	1299	WRITE(knum,*) 'Number of subdomains located along the north fold : ', ndim_rank_north
	1300	WRITE(knum,*) 'Rank of the subdomains located along the north fold : ', ndim_rank_north
	1301	DO jn = 1, ndim_rank_north, 5
	1302	WRITE(knum,*) nrank_north( jn:MINVAL( (/jn+4,ndim_rank_north/) ) )
	1303	END DO
	1304	ENDIF
	1305
	1306	nfd_nbnei = 0 ! defaul def (useless?)
	1307	IF( ln_nnogather ) THEN
[9436]	1308	!
[15267]	1309	! Use the "gather nfd" to know how to do the nfd: for ji point, which process send data from which of its ji-index?
	1310	! Note that nfd is perfectly symetric: I receive data from X <=> I send data to X (-> no deadlock)
[9436]	1311	!
[15267]	1312	zinfo(:,:,:,:) = HUGE(0._wp) ! default def to make sure we don't use the halos
	1313	DO jg = 1, 4 ! grid type: T, U, V, F
	1314	DO jj = nn_hls+1, jpj-nn_hls ! inner domain (warning do_loop_substitute not yet defined)
	1315	DO ji = nn_hls+1, jpi-nn_hls ! inner domain (warning do_loop_substitute not yet defined)
	1316	zinfo(ji,jj,1,jg) = REAL(narea, wp) ! mpi_rank + 1 (as default lbc_lnk fill is 0
	1317	zinfo(ji,jj,2,jg) = REAL(ji, wp) ! ji of this proc
	1318	END DO
	1319	END DO
	1320	END DO
[9436]	1321	!
[15267]	1322	ln_nnogather = .FALSE. ! force "classical" North pole folding to fill all halos -> should be no more HUGE values...
	1323	CALL lbc_lnk( 'mppini', zinfo(:,:,:,1), 'T', 1._wp ) ! Do 4 calls instead of 1 to save memory as the nogather version
	1324	CALL lbc_lnk( 'mppini', zinfo(:,:,:,2), 'U', 1._wp ) ! creates buffer arrays with jpiglo as the first dimension
	1325	CALL lbc_lnk( 'mppini', zinfo(:,:,:,3), 'V', 1._wp ) !
	1326	CALL lbc_lnk( 'mppini', zinfo(:,:,:,4), 'F', 1._wp ) !
	1327	ln_nnogather = .TRUE.
	1328
	1329	IF( l_IdoNFold ) THEN ! only the procs involed in the NFD must take care of this
	1330
	1331	ALLOCATE( nfd_rksnd(jpi,4), nfd_jisnd(jpi,4) ) ! neighbour rand and remote ji-index for each grid (T, U, V, F)
	1332	nfd_rksnd(:,:) = NINT( zinfo(:, jpj, 1, :) ) - 1 ! neighbour MPI rank
	1333	nfd_jisnd(:,:) = NINT( zinfo(:, jpj, 2, :) ) - nn_hls ! neighbour ji index (shifted as we don't send the halos)
	1334	WHERE( nfd_rksnd == -1 ) nfd_jisnd = 1 ! use ji=1 if no neighbour, see mpp_nfd_generic.h90
	1335
	1336	nfd_nbnei = 1 ! Number of neighbour sending data for the nfd. We have at least 1 neighbour!
	1337	irknei(1) = nfd_rksnd(1,1) ! which is the 1st one (I can be neighbour of myself, exclude land-proc are also ok)
	1338	DO jg = 1, 4
	1339	DO ji = 1, jpi ! we must be able to fill the full line including halos
	1340	lnew = .TRUE. ! new neighbour?
	1341	DO jn = 1, nfd_nbnei
	1342	IF( irknei(jn) == nfd_rksnd(ji,jg) ) lnew = .FALSE. ! already found
	1343	END DO
	1344	IF( lnew ) THEN
	1345	jn = nfd_nbnei + 1
	1346	nfd_nbnei = jn
	1347	irknei(jn) = nfd_rksnd(ji,jg)
	1348	ENDIF
	1349	END DO
	1350	END DO
	1351
	1352	ALLOCATE( nfd_rknei(nfd_nbnei) )
	1353	nfd_rknei(:) = irknei(1:nfd_nbnei)
	1354	! re-number nfd_rksnd according to the indexes of nfd_rknei
	1355	DO jn = 1, nfd_nbnei
	1356	WHERE( nfd_rksnd == nfd_rknei(jn) ) nfd_rksnd = jn
	1357	END DO
	1358
	1359	IF( ldwrtlay ) THEN
	1360	WRITE(knum,*)
	1361	WRITE(knum,*) 'north fold exchanges with explicit point-to-point messaging :'
	1362	WRITE(knum,*) ' number of neighbours for the NF: nfd_nbnei : ', nfd_nbnei
	1363	IF( nfd_nbnei > 0 ) WRITE(knum,*) ' neighbours MPI ranks : ', nfd_rknei
	1364	ENDIF
	1365
	1366	ENDIF ! l_IdoNFold
	1367	!
	1368	ENDIF ! ln_nnogather
[9436]	1369	!
[13286]	1370	END SUBROUTINE init_nfdcom
[9436]	1371
[88]	1372
[13286]	1373	SUBROUTINE init_doloop
	1374	!!----------------------------------------------------------------------
	1375	!! * ROUTINE init_doloop *
	1376	!!
	1377	!! ** Purpose : set the starting/ending indices of DO-loop
	1378	!! These indices are used in do_loop_substitute.h90
	1379	!!----------------------------------------------------------------------
	1380	!
[14433]	1381	Nis0 = 1+nn_hls
	1382	Njs0 = 1+nn_hls
	1383	Nie0 = jpi-nn_hls
	1384	Nje0 = jpj-nn_hls
[14072]	1385	!
[13286]	1386	Ni_0 = Nie0 - Nis0 + 1
	1387	Nj_0 = Nje0 - Njs0 + 1
	1388	!
[14433]	1389	jpkm1 = jpk-1 ! " "
	1390	!
[13286]	1391	END SUBROUTINE init_doloop
[14072]	1392
[15267]	1393
	1394	SUBROUTINE init_locglo
	1395	!!----------------------------------------------------------------------
	1396	!! * ROUTINE init_locglo *
	1397	!!
	1398	!! ** Purpose : initialization of global domain <--> local domain indices
	1399	!!
	1400	!! ** Method :
	1401	!!
	1402	!! ** Action : - mig , mjg : local domain indices ==> global domain, including halos, indices
	1403	!! - mig0, mjg0: local domain indices ==> global domain, excluding halos, indices
	1404	!! - mi0 , mi1 : global domain indices ==> local domain indices
	1405	!! - mj0 , mj1 (if global point not in the local domain ==> mi0>mi1 and/or mj0>mj1)
	1406	!!----------------------------------------------------------------------
	1407	INTEGER :: ji, jj ! dummy loop argument
	1408	!!----------------------------------------------------------------------
	1409	!
	1410	ALLOCATE( mig(jpi), mjg(jpj), mig0(jpi), mjg0(jpj) )
	1411	ALLOCATE( mi0(jpiglo), mi1(jpiglo), mj0(jpjglo), mj1(jpjglo) )
	1412	!
	1413	DO ji = 1, jpi ! local domain indices ==> global domain indices, including halos
	1414	mig(ji) = ji + nimpp - 1
	1415	END DO
	1416	DO jj = 1, jpj
	1417	mjg(jj) = jj + njmpp - 1
	1418	END DO
	1419	! ! local domain indices ==> global domain indices, excluding halos
	1420	!
	1421	mig0(:) = mig(:) - nn_hls
	1422	mjg0(:) = mjg(:) - nn_hls
	1423	! ! global domain, including halos, indices ==> local domain indices
	1424	! ! (return (m.0,m.1)=(1,0) if data domain gridpoint is to the west/south of the
	1425	! ! local domain, or (m.0,m.1)=(jp.+1,jp.) to the east/north of local domain.
	1426	DO ji = 1, jpiglo
	1427	mi0(ji) = MAX( 1 , MIN( ji - nimpp + 1, jpi+1 ) )
	1428	mi1(ji) = MAX( 0 , MIN( ji - nimpp + 1, jpi ) )
	1429	END DO
	1430	DO jj = 1, jpjglo
	1431	mj0(jj) = MAX( 1 , MIN( jj - njmpp + 1, jpj+1 ) )
	1432	mj1(jj) = MAX( 0 , MIN( jj - njmpp + 1, jpj ) )
	1433	END DO
	1434	!
	1435	END SUBROUTINE init_locglo
	1436
[3]	1437	!!======================================================================
	1438	END MODULE mppini

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

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

Context Navigation

source: NEMO/trunk/src/OCE/LBC/mppini.F90

Download in other formats: