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/trunk/src/OCE/LBC – NEMO

Context Navigation

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

Last change on this file since 15033 was 15033, checked in by smasson, 3 years ago
trunk: suppress jpim1 et jpjm1, #2699
Property svn:keywords set to `Id`
File size: 65.7 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	!
[14072]	25	USE lbcnfd , ONLY : isendto, nsndto ! 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'
	170	IF( jpni < 1 .OR. jpnj < 1 ) THEN
	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	!
[14433]	333	CALL init_doloop ! set start/end indices of do-loop, depending on the halo width value (nn_hls)
	334	!
[9019]	335	IF(lwp) THEN
	336	WRITE(numout,*)
[10425]	337	WRITE(numout,*) 'MPI Message Passing MPI - domain lay out over processors'
	338	WRITE(numout,*)
[9019]	339	WRITE(numout,*) ' defines mpp subdomains'
[14072]	340	WRITE(numout,*) ' jpni = ', jpni
[10425]	341	WRITE(numout,*) ' jpnj = ', jpnj
[13286]	342	WRITE(numout,*) ' jpnij = ', jpnij
[14433]	343	WRITE(numout,*) ' nimpp = ', nimpp
	344	WRITE(numout,*) ' njmpp = ', njmpp
[9019]	345	WRITE(numout,*)
[13286]	346	WRITE(numout,*) ' sum ijpi(i,1) = ', sum(ijpi(:,1)), ' jpiglo = ', jpiglo
[14433]	347	WRITE(numout,*) ' sum ijpj(1,j) = ', SUM(ijpj(1,:)), ' jpjglo = ', jpjglo
	348
	349	! Subdomain grid print
[10425]	350	ifreq = 4
	351	il1 = 1
	352	DO jn = 1, (jpni-1)/ifreq+1
	353	il2 = MIN(jpni,il1+ifreq-1)
	354	WRITE(numout,*)
	355	WRITE(numout,9400) ('***',ji=il1,il2-1)
	356	DO jj = jpnj, 1, -1
	357	WRITE(numout,9403) (' ',ji=il1,il2-1)
[13286]	358	WRITE(numout,9402) jj, (ijpi(ji,jj),ijpj(ji,jj),ji=il1,il2)
[10425]	359	WRITE(numout,9404) (ipproc(ji,jj),ji=il1,il2)
	360	WRITE(numout,9403) (' ',ji=il1,il2-1)
	361	WRITE(numout,9400) ('***',ji=il1,il2-1)
	362	END DO
	363	WRITE(numout,9401) (ji,ji=il1,il2)
	364	il1 = il1+ifreq
	365	END DO
	366	9400 FORMAT(' *' ,20('***********',a3) )
	367	9403 FORMAT(' * ',20(' * ',a3) )
	368	9401 FORMAT(' ' ,20(' ',i3,' ') )
	369	9402 FORMAT(' ',i3,' * ',20(i3,' x',i3,' * ') )
[11640]	370	9404 FORMAT(' * ' ,20(' ' ,i4,' * ') )
[10425]	371	ENDIF
[14433]	372	!
	373	! Store informations for the north pole folding communications
	374	nfproc(:) = ipproc(:,jpnj)
	375	nfimpp(:) = iimppt(:,jpnj)
	376	nfjpi (:) = ijpi(:,jpnj)
	377	!
	378	! 3. Define Western, Eastern, Southern and Northern neighbors + corners in the subdomain grid reference
	379	! ------------------------------------------------------------------------------------------------------
	380	!
	381	! note that North fold is has specific treatment for its MPI communications.
	382	! This must not be treated as a "usual" communication with a northern neighbor.
	383	! -> North fold processes have no Northern neighbor in the definition done bellow
	384	!
	385	llmpi_Iperio = jpni > 1 .AND. l_Iperio ! do i-periodicity with an MPI communication?
	386	llmpi_Jperio = jpnj > 1 .AND. l_Jperio ! do j-periodicity with an MPI communication?
	387	!
	388	l_SelfPerio(1:2) = l_Iperio .AND. jpni == 1 ! west, east periodicity by itself
	389	l_SelfPerio(3:4) = l_Jperio .AND. jpnj == 1 ! south, north periodicity by itself
	390	l_SelfPerio(5:8) = l_SelfPerio(jpwe) .AND. l_SelfPerio(jpso) ! corners bi-periodicity by itself
	391	!
	392	! define neighbors mapping (1/2): default definition: ignore if neighbours are land-only subdomains or not
	393	DO jj = 1, jpnj
	394	DO ji = 1, jpni
	395	!
	396	IF ( llisOce(ji,jj) ) THEN ! this subdomain has some ocean: it has neighbours
	397	!
	398	inum0 = ji - 1 + ( jj - 1 ) * jpni ! index in the subdomains grid. start at 0
	399	!
	400	! Is there a neighbor?
	401	llnei(jpwe,ji,jj) = ji > 1 .OR. llmpi_Iperio ! West nei exists if not the first column or llmpi_Iperio
	402	llnei(jpea,ji,jj) = ji < jpni .OR. llmpi_Iperio ! East nei exists if not the last column or llmpi_Iperio
	403	llnei(jpso,ji,jj) = jj > 1 .OR. llmpi_Jperio ! South nei exists if not the first line or llmpi_Jperio
	404	llnei(jpno,ji,jj) = jj < jpnj .OR. llmpi_Jperio ! North nei exists if not the last line or llmpi_Jperio
	405	llnei(jpsw,ji,jj) = llnei(jpwe,ji,jj) .AND. llnei(jpso,ji,jj) ! So-We nei exists if both South and West nei exist
	406	llnei(jpse,ji,jj) = llnei(jpea,ji,jj) .AND. llnei(jpso,ji,jj) ! So-Ea nei exists if both South and East nei exist
	407	llnei(jpnw,ji,jj) = llnei(jpwe,ji,jj) .AND. llnei(jpno,ji,jj) ! No-We nei exists if both North and West nei exist
	408	llnei(jpne,ji,jj) = llnei(jpea,ji,jj) .AND. llnei(jpno,ji,jj) ! No-Ea nei exists if both North and East nei exist
	409	!
	410	! Which index (starting at 0) have neighbors in the subdomains grid?
	411	IF( llnei(jpwe,ji,jj) ) inei(jpwe,ji,jj) = inum0 - 1 + jpni * COUNT( (/ ji == 1 /) )
	412	IF( llnei(jpea,ji,jj) ) inei(jpea,ji,jj) = inum0 + 1 - jpni * COUNT( (/ ji == jpni /) )
	413	IF( llnei(jpso,ji,jj) ) inei(jpso,ji,jj) = inum0 - jpni + jpni * jpnj * COUNT( (/ jj == 1 /) )
	414	IF( llnei(jpno,ji,jj) ) inei(jpno,ji,jj) = inum0 + jpni - jpni * jpnj * COUNT( (/ jj == jpnj /) )
	415	IF( llnei(jpsw,ji,jj) ) inei(jpsw,ji,jj) = inei(jpso,ji,jj) - 1 + jpni * COUNT( (/ ji == 1 /) )
	416	IF( llnei(jpse,ji,jj) ) inei(jpse,ji,jj) = inei(jpso,ji,jj) + 1 - jpni * COUNT( (/ ji == jpni /) )
	417	IF( llnei(jpnw,ji,jj) ) inei(jpnw,ji,jj) = inei(jpno,ji,jj) - 1 + jpni * COUNT( (/ ji == 1 /) )
	418	IF( llnei(jpne,ji,jj) ) inei(jpne,ji,jj) = inei(jpno,ji,jj) + 1 - jpni * COUNT( (/ ji == jpni /) )
	419	!
	420	ELSE ! land-only domain has no neighbour
	421	llnei(:,ji,jj) = .FALSE.
	422	ENDIF
	423	!
	424	END DO
[9019]	425	END DO
	426	!
[14433]	427	! define neighbors mapping (2/2): check if neighbours are not land-only subdomains
	428	DO jj = 1, jpnj
	429	DO ji = 1, jpni
	430	DO jn = 1, 8
	431	IF( llnei(jn,ji,jj) ) THEN ! if a neighbour is existing -> this should not be a land-only domain
	432	ii = 1 + MOD( inei(jn,ji,jj) , jpni )
	433	ij = 1 + inei(jn,ji,jj) / jpni
	434	llnei(jn,ji,jj) = llisOce( ii, ij )
	435	ENDIF
	436	END DO
	437	END DO
	438	END DO
[13286]	439	!
[14433]	440	! update index of the neighbours in the subdomains grid
	441	WHERE( .NOT. llnei ) inei = -1
[13286]	442	!
[9019]	443	! Save processor layout in ascii file
[10570]	444	IF (llwrtlay) THEN
[7646]	445	CALL ctl_opn( inum, 'layout.dat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE., narea )
[14433]	446	WRITE(inum,'(a)') ' jpnij jpimax jpjmax jpk jpiglo jpjglo ( local: narea jpi jpj )'
	447	WRITE(inum,'(6i7,a,3i7,a)') jpnij,jpimax,jpjmax,jpk,jpiglo,jpjglo,' ( local: ',narea,jpi,jpj,' )'
	448	WRITE(inum,*)
	449	WRITE(inum, *) '------------------------------------'
	450	WRITE(inum,'(a,i2)') ' Mapping of the default neighnourgs '
	451	WRITE(inum, *) '------------------------------------'
	452	WRITE(inum,*)
	453	WRITE(inum,'(a)') ' rank ii ij jpi jpj nimpp njmpp mpiwe mpiea mpiso mpino mpisw mpise mpinw mpine'
	454	DO jp = 1, jpnij
	455	ii = iin(jp)
	456	ij = ijn(jp)
	457	WRITE(inum,'(15i6)') jp-1, ii, ij, ijpi(ii,ij), ijpj(ii,ij), iimppt(ii,ij), ijmppt(ii,ij), inei(:,ii,ij)
	458	END DO
	459	ENDIF
[9019]	460
[14433]	461	!
	462	! 4. Define Western, Eastern, Southern and Northern neighbors + corners for each mpi process
	463	! ------------------------------------------------------------------------------------------
	464	!
	465	! rewrite information from "subdomain grid" to mpi process list
	466	! Warning, for example:
	467	! position of the northern neighbor in the "subdomain grid"
	468	! position of the northern neighbor in the "mpi process list"
	469
	470	! default definition: no neighbors
	471	impi(:,:) = -1 ! (starting at 0, -1 if no neighbourg)
	472
	473	DO jp = 1, jpnij
	474	ii = iin(jp)
	475	ij = ijn(jp)
	476	DO jn = 1, 8
	477	IF( llnei(jn,ii,ij) ) THEN ! must be tested as some land-domain can be kept to fit mppsize
	478	ii2 = 1 + MOD( inei(jn,ii,ij) , jpni )
	479	ij2 = 1 + inei(jn,ii,ij) / jpni
	480	impi(jn,jp) = ipproc( ii2, ij2 )
	481	ENDIF
[7646]	482	END DO
[14433]	483	END DO
[3]	484
[9019]	485	!
[14433]	486	! 4. keep information for the local process
	487	! -----------------------------------------
	488	!
	489	! set default neighbours
	490	mpinei(:) = impi(:,narea)
	491	DO jh = 1, n_hlsmax
	492	mpiSnei(jh,:) = impi(:,narea) ! default definition
	493	mpiRnei(jh,:) = impi(:,narea)
	494	END DO
	495	!
[3]	496	IF(lwp) THEN
[6412]	497	WRITE(numout,*)
[9169]	498	WRITE(numout,*) ' resulting internal parameters : '
[14433]	499	WRITE(numout,*) ' narea = ', narea
	500	WRITE(numout,*) ' mpi nei west = ', mpinei(jpwe) , ' mpi nei east = ', mpinei(jpea)
	501	WRITE(numout,*) ' mpi nei south = ', mpinei(jpso) , ' mpi nei north = ', mpinei(jpno)
	502	WRITE(numout,*) ' mpi nei so-we = ', mpinei(jpsw) , ' mpi nei so-ea = ', mpinei(jpse)
	503	WRITE(numout,*) ' mpi nei no-we = ', mpinei(jpnw) , ' mpi nei no-ea = ', mpinei(jpne)
[3]	504	ENDIF
[9019]	505	! ! Prepare mpp north fold
[14433]	506	!
	507	llmpiNFold = jpni > 1 .AND. l_NFold ! is the North fold done with an MPI communication?
	508	l_IdoNFold = ijn(narea) == jpnj .AND. l_NFold ! is this process doing North fold?
	509	!
	510	IF( llmpiNFold ) THEN
[3]	511	CALL mpp_ini_north
[9917]	512	IF (lwp) THEN
	513	WRITE(numout,*)
	514	WRITE(numout,*) ' ==>>> North fold boundary prepared for jpni >1'
[10570]	515	ENDIF
[14433]	516	IF (llwrtlay) THEN ! additional prints in layout.dat
[9917]	517	WRITE(inum,*)
	518	WRITE(inum,*)
[14433]	519	WRITE(inum,*) 'Number of subdomains located along the north fold : ', ndim_rank_north
[9917]	520	WRITE(inum,*) 'Rank of the subdomains located along the north fold : ', ndim_rank_north
[14433]	521	DO jp = 1, ndim_rank_north, 5
	522	WRITE(inum,*) nrank_north( jp:MINVAL( (/jp+4,ndim_rank_north/) ) )
[9917]	523	END DO
	524	ENDIF
[14433]	525	IF ( l_IdoNFold .AND. ln_nnogather ) THEN
	526	CALL init_nfdcom ! northfold neighbour lists
	527	IF (llwrtlay) THEN
	528	WRITE(inum,*)
	529	WRITE(inum,*) 'north fold exchanges with explicit point-to-point messaging :'
	530	WRITE(inum,*) ' nsndto : ', nsndto
	531	WRITE(inum,*) ' isendto : ', isendto(1:nsndto)
	532	ENDIF
	533	ENDIF
[6412]	534	ENDIF
[9019]	535	!
[14433]	536	CALL mpp_ini_nc(nn_hls) ! Initialize communicator for neighbourhood collective communications
	537	DO jh = 1, n_hlsmax
	538	mpi_nc_com4(jh) = mpi_nc_com4(nn_hls) ! default definition
	539	mpi_nc_com8(jh) = mpi_nc_com8(nn_hls)
	540	END DO
[13982]	541	!
[15023]	542	IF( jpnij > 1) CALL init_excl_landpt ! exclude exchanges which contain only land points
[14072]	543	!
[14433]	544	! Save processor layout changes in ascii file
	545	DO jh = 1, n_hlsmax ! different halo size
	546	DO ji = 1, 8
	547	ichanged(16*(jh-1) +ji) = COUNT( mpinei(ji:ji) /= mpiSnei(jh,ji:ji) )
	548	ichanged(16*(jh-1)+8+ji) = COUNT( mpinei(ji:ji) /= mpiRnei(jh,ji:ji) )
	549	END DO
	550	END DO
	551	CALL mpp_sum( "mpp_init", ichanged ) ! must be called by all processes
	552	IF (llwrtlay) THEN
	553	WRITE(inum,*)
	554	WRITE(inum, *) '----------------------------------------------------------------------'
	555	WRITE(inum,'(a,i2)') ' Mapping of the neighnourgs once excluding comm with only land points '
	556	WRITE(inum, *) '----------------------------------------------------------------------'
	557	DO jh = 1, n_hlsmax ! different halo size
	558	WRITE(inum,*)
	559	WRITE(inum,'(a,i2)') 'halo size: ', jh
	560	WRITE(inum, *) '---------'
	561	WRITE(inum,'(a)') ' rank ii ij mpiwe mpiea mpiso mpino mpisw mpise mpinw mpine'
	562	WRITE(inum, '(11i6,a)') narea-1, iin(narea), ijn(narea), mpinei(:), ' <- Org'
	563	WRITE(inum,'(18x,8i6,a,i1,a)') mpiSnei(jh,:), ' <- Send ', COUNT( mpinei(:) /= mpiSnei(jh,:) ), ' modif'
	564	WRITE(inum,'(18x,8i6,a,i1,a)') mpiRnei(jh,:), ' <- Recv ', COUNT( mpinei(:) /= mpiRnei(jh,:) ), ' modif'
	565	WRITE(inum,*) ' total changes among all mpi tasks:'
	566	WRITE(inum,*) ' mpiwe mpiea mpiso mpino mpisw mpise mpinw mpine'
	567	WRITE(inum,'(a,8i6)') ' Send: ', ichanged(jh16-15:jh16-8)
	568	WRITE(inum,'(a,8i6)') ' Recv: ', ichanged(jh16 -7:jh16 )
	569	END DO
[9917]	570	ENDIF
[9436]	571	!
[14433]	572	CALL init_ioipsl ! Prepare NetCDF output file (if necessary)
	573	!
[14072]	574	IF (llwrtlay) CLOSE(inum)
[9917]	575	!
[14433]	576	DEALLOCATE(iin, ijn, iimppt, ijmppt, ijpi, ijpj, ipproc, inei, llnei, impi, llisOce)
[9444]	577	!
[9019]	578	END SUBROUTINE mpp_init
[3]	579
[13438]	580	#endif
[3]	581
[13286]	582	SUBROUTINE mpp_basesplit( kiglo, kjglo, khls, knbi, knbj, kimax, kjmax, kimppt, kjmppt, klci, klcj)
[9019]	583	!!----------------------------------------------------------------------
[13286]	584	!! * ROUTINE mpp_basesplit *
[14072]	585	!!
[10425]	586	!! ** Purpose : Lay out the global domain over processors.
[9019]	587	!!
[10425]	588	!! ** Method : Global domain is distributed in smaller local domains.
[9019]	589	!!
[10425]	590	!! ** Action : - set for all knbi*knbj domains:
	591	!! kimppt : longitudinal index
	592	!! kjmppt : latitudinal index
	593	!! klci : first dimension
	594	!! klcj : second dimension
[9019]	595	!!----------------------------------------------------------------------
[13286]	596	INTEGER, INTENT(in ) :: kiglo, kjglo
	597	INTEGER, INTENT(in ) :: khls
[10425]	598	INTEGER, INTENT(in ) :: knbi, knbj
	599	INTEGER, INTENT( out) :: kimax, kjmax
	600	INTEGER, DIMENSION(knbi,knbj), OPTIONAL, INTENT( out) :: kimppt, kjmppt
	601	INTEGER, DIMENSION(knbi,knbj), OPTIONAL, INTENT( out) :: klci, klcj
	602	!
	603	INTEGER :: ji, jj
[14072]	604	INTEGER :: i2hls
[10542]	605	INTEGER :: iresti, irestj, irm, ijpjmin
[9019]	606	!!----------------------------------------------------------------------
[13286]	607	i2hls = 2*khls
[9019]	608	!
[10425]	609	#if defined key_nemocice_decomp
[13286]	610	kimax = ( nx_global+2-i2hls + (knbi-1) ) / knbi + i2hls ! first dim.
[14072]	611	kjmax = ( ny_global+2-i2hls + (knbj-1) ) / knbj + i2hls ! second dim.
[10425]	612	#else
[13286]	613	kimax = ( kiglo - i2hls + (knbi-1) ) / knbi + i2hls ! first dim.
	614	kjmax = ( kjglo - i2hls + (knbj-1) ) / knbj + i2hls ! second dim.
[10425]	615	#endif
	616	IF( .NOT. PRESENT(kimppt) ) RETURN
[9019]	617	!
[10425]	618	! 1. Dimension arrays for subdomains
	619	! -----------------------------------
	620	! Computation of local domain sizes klci() klcj()
[13286]	621	! These dimensions depend on global sizes knbi,knbj and kiglo,kjglo
[10425]	622	! The subdomains are squares lesser than or equal to the global
	623	! dimensions divided by the number of processors minus the overlap array.
[9019]	624	!
[13286]	625	iresti = 1 + MOD( kiglo - i2hls - 1 , knbi )
	626	irestj = 1 + MOD( kjglo - i2hls - 1 , knbj )
[9168]	627	!
[10425]	628	! Need to use kimax and kjmax here since jpi and jpj not yet defined
	629	#if defined key_nemocice_decomp
	630	! Change padding to be consistent with CICE
[13286]	631	klci(1:knbi-1,: ) = kimax
	632	klci( knbi ,: ) = kiglo - (knbi - 1) * (kimax - i2hls)
	633	klcj(: ,1:knbj-1) = kjmax
	634	klcj(: , knbj ) = kjglo - (knbj - 1) * (kjmax - i2hls)
[10425]	635	#else
	636	klci(1:iresti ,:) = kimax
	637	klci(iresti+1:knbi ,:) = kimax-1
[14848]	638	IF( MINVAL(klci) < 3*khls ) THEN
	639	WRITE(ctmp1,) ' mpp_basesplit: minimum value of jpi must be >= ', 3khls
[10542]	640	WRITE(ctmp2,*) ' We have ', MINVAL(klci)
[14848]	641	CALL ctl_stop( 'STOP', ctmp1, ctmp2 )
[10542]	642	ENDIF
[14433]	643	IF( l_NFold ) THEN
[10542]	644	! minimize the size of the last row to compensate for the north pole folding coast
[14433]	645	IF( c_NFtype == 'T' ) ijpjmin = 2+3*khls ! V and F folding must be outside of southern halos
	646	IF( c_NFtype == 'F' ) ijpjmin = 1+3*khls ! V and F folding must be outside of southern halos
	647	irm = knbj - irestj ! total number of lines to be removed
	648	klcj(:,knbj) = MAX( ijpjmin, kjmax-irm ) ! we must have jpj >= ijpjmin in the last row
	649	irm = irm - ( kjmax - klcj(1,knbj) ) ! remaining number of lines to remove
[13286]	650	irestj = knbj - 1 - irm
[10542]	651	klcj(:, irestj+1:knbj-1) = kjmax-1
	652	ELSE
[13286]	653	klcj(:, irestj+1:knbj ) = kjmax-1
[10542]	654	ENDIF
[13286]	655	klcj(:,1:irestj) = kjmax
[14848]	656	IF( MINVAL(klcj) < 3*khls ) THEN
	657	WRITE(ctmp1,) ' mpp_basesplit: minimum value of jpj must be >= ', 3khls
[10542]	658	WRITE(ctmp2,*) ' We have ', MINVAL(klcj)
	659	CALL ctl_stop( 'STOP', ctmp1, ctmp2 )
	660	ENDIF
[10425]	661	#endif
[3]	662
[10425]	663	! 2. Index arrays for subdomains
	664	! -------------------------------
	665	kimppt(:,:) = 1
	666	kjmppt(:,:) = 1
	667	!
	668	IF( knbi > 1 ) THEN
	669	DO jj = 1, knbj
	670	DO ji = 2, knbi
[13286]	671	kimppt(ji,jj) = kimppt(ji-1,jj) + klci(ji-1,jj) - i2hls
[10425]	672	END DO
	673	END DO
[9019]	674	ENDIF
	675	!
[10425]	676	IF( knbj > 1 )THEN
	677	DO jj = 2, knbj
	678	DO ji = 1, knbi
[13286]	679	kjmppt(ji,jj) = kjmppt(ji,jj-1) + klcj(ji,jj-1) - i2hls
[10425]	680	END DO
	681	END DO
	682	ENDIF
[14072]	683
[13286]	684	END SUBROUTINE mpp_basesplit
[9019]	685
	686
[13286]	687	SUBROUTINE bestpartition( knbij, knbi, knbj, knbcnt, ldlist )
[10425]	688	!!----------------------------------------------------------------------
[13286]	689	!! * ROUTINE bestpartition *
[10425]	690	!!
	691	!! ** Purpose :
	692	!!
	693	!! ** Method :
	694	!!----------------------------------------------------------------------
[13490]	695	INTEGER, INTENT(in ) :: knbij ! total number of subdomains (knbi*knbj)
[10425]	696	INTEGER, OPTIONAL, INTENT( out) :: knbi, knbj ! number if subdomains along i and j (knbi and knbj)
	697	INTEGER, OPTIONAL, INTENT( out) :: knbcnt ! number of land subdomains
	698	LOGICAL, OPTIONAL, INTENT(in ) :: ldlist ! .true.: print the list the best domain decompositions (with land)
	699	!
	700	INTEGER :: ji, jj, ii, iitarget
	701	INTEGER :: iszitst, iszjtst
	702	INTEGER :: isziref, iszjref
[13490]	703	INTEGER :: iszimin, iszjmin
[10425]	704	INTEGER :: inbij, iszij
[11536]	705	INTEGER :: inbimax, inbjmax, inbijmax, inbijold
[10425]	706	INTEGER :: isz0, isz1
	707	INTEGER, DIMENSION( :), ALLOCATABLE :: indexok
	708	INTEGER, DIMENSION( :), ALLOCATABLE :: inbi0, inbj0, inbij0 ! number of subdomains along i,j
	709	INTEGER, DIMENSION( :), ALLOCATABLE :: iszi0, iszj0, iszij0 ! max size of the subdomains along i,j
	710	INTEGER, DIMENSION( :), ALLOCATABLE :: inbi1, inbj1, inbij1 ! number of subdomains along i,j
	711	INTEGER, DIMENSION( :), ALLOCATABLE :: iszi1, iszj1, iszij1 ! max size of the subdomains along i,j
	712	LOGICAL :: llist
	713	LOGICAL, DIMENSION(:,:), ALLOCATABLE :: llmsk2d ! max size of the subdomains along i,j
[14433]	714	LOGICAL, DIMENSION(:,:), ALLOCATABLE :: llisOce ! - -
[10425]	715	REAL(wp):: zpropland
	716	!!----------------------------------------------------------------------
	717	!
	718	llist = .FALSE.
	719	IF( PRESENT(ldlist) ) llist = ldlist
	720
	721	CALL mpp_init_landprop( zpropland ) ! get the proportion of land point over the gloal domain
	722	inbij = NINT( REAL(knbij, wp) / ( 1.0 - zpropland ) ) ! define the largest possible value for jpni*jpnj
	723	!
	724	IF( llist ) THEN ; inbijmax = inbij*2
	725	ELSE ; inbijmax = inbij
	726	ENDIF
	727	!
	728	ALLOCATE(inbi0(inbijmax),inbj0(inbijmax),iszi0(inbijmax),iszj0(inbijmax))
	729	!
	730	inbimax = 0
	731	inbjmax = 0
[13490]	732	isziref = jpiglo*jpjglo+1 ! define a value that is larger than the largest possible
	733	iszjref = jpiglo*jpjglo+1
[10425]	734	!
[14848]	735	iszimin = 3*nn_hls ! minimum size of the MPI subdomain so halos are always adressing neighbor inner domain
	736	iszjmin = 3*nn_hls
[14433]	737	IF( c_NFtype == 'T' ) iszjmin = MAX(iszjmin, 2+3*nn_hls) ! V and F folding must be outside of southern halos
	738	IF( c_NFtype == 'F' ) iszjmin = MAX(iszjmin, 1+3*nn_hls) ! V and F folding must be outside of southern halos
[13490]	739	!
[10425]	740	! get the list of knbi that gives a smaller jpimax than knbi-1
	741	! get the list of knbj that gives a smaller jpjmax than knbj-1
[14072]	742	DO ji = 1, inbijmax
[10425]	743	#if defined key_nemocice_decomp
	744	iszitst = ( nx_global+2-2nn_hls + (ji-1) ) / ji + 2nn_hls ! first dim.
	745	#else
[13490]	746	iszitst = ( Ni0glo + (ji-1) ) / ji + 2*nn_hls ! max subdomain i-size
[10425]	747	#endif
[13490]	748	IF( iszitst < isziref .AND. iszitst >= iszimin ) THEN
[10425]	749	isziref = iszitst
	750	inbimax = inbimax + 1
	751	inbi0(inbimax) = ji
	752	iszi0(inbimax) = isziref
	753	ENDIF
	754	#if defined key_nemocice_decomp
	755	iszjtst = ( ny_global+2-2nn_hls + (ji-1) ) / ji + 2nn_hls ! first dim.
	756	#else
[13490]	757	iszjtst = ( Nj0glo + (ji-1) ) / ji + 2*nn_hls ! max subdomain j-size
[10425]	758	#endif
[13490]	759	IF( iszjtst < iszjref .AND. iszjtst >= iszjmin ) THEN
[10425]	760	iszjref = iszjtst
	761	inbjmax = inbjmax + 1
	762	inbj0(inbjmax) = ji
	763	iszj0(inbjmax) = iszjref
	764	ENDIF
	765	END DO
[14848]	766	IF( inbimax == 0 ) THEN
	767	WRITE(ctmp1,'(a,i2,a,i2)') ' mpp_ini bestpartition: Ni0glo (', Ni0glo, ') is too small to be used with nn_hls = ', nn_hls
	768	CALL ctl_stop( 'STOP', ctmp1 )
	769	ENDIF
	770	IF( inbjmax == 0 ) THEN
	771	WRITE(ctmp1,'(a,i2,a,i2)') ' mpp_ini bestpartition: Nj0glo (', Nj0glo, ') is too small to be used with nn_hls = ', nn_hls
	772	CALL ctl_stop( 'STOP', ctmp1 )
	773	ENDIF
[10425]	774
	775	! combine these 2 lists to get all possible knbi*knbj < inbijmax
	776	ALLOCATE( llmsk2d(inbimax,inbjmax) )
	777	DO jj = 1, inbjmax
	778	DO ji = 1, inbimax
	779	IF ( inbi0(ji) * inbj0(jj) <= inbijmax ) THEN ; llmsk2d(ji,jj) = .TRUE.
	780	ELSE ; llmsk2d(ji,jj) = .FALSE.
	781	ENDIF
	782	END DO
	783	END DO
	784	isz1 = COUNT(llmsk2d)
	785	ALLOCATE( inbi1(isz1), inbj1(isz1), iszi1(isz1), iszj1(isz1) )
	786	ii = 0
	787	DO jj = 1, inbjmax
	788	DO ji = 1, inbimax
	789	IF( llmsk2d(ji,jj) .EQV. .TRUE. ) THEN
	790	ii = ii + 1
	791	inbi1(ii) = inbi0(ji)
	792	inbj1(ii) = inbj0(jj)
	793	iszi1(ii) = iszi0(ji)
	794	iszj1(ii) = iszj0(jj)
[14433]	795	ENDIF
[10425]	796	END DO
	797	END DO
	798	DEALLOCATE( inbi0, inbj0, iszi0, iszj0 )
	799	DEALLOCATE( llmsk2d )
	800
	801	ALLOCATE( inbij1(isz1), iszij1(isz1) )
	802	inbij1(:) = inbi1(:) * inbj1(:)
	803	iszij1(:) = iszi1(:) * iszj1(:)
	804
[13286]	805	! if there is no land and no print
[10425]	806	IF( .NOT. llist .AND. numbot == -1 .AND. numbdy == -1 ) THEN
	807	! get the smaller partition which gives the smallest subdomain size
	808	ii = MINLOC(inbij1, mask = iszij1 == MINVAL(iszij1), dim = 1)
	809	knbi = inbi1(ii)
	810	knbj = inbj1(ii)
	811	IF(PRESENT(knbcnt)) knbcnt = 0
	812	DEALLOCATE( inbi1, inbj1, inbij1, iszi1, iszj1, iszij1 )
	813	RETURN
	814	ENDIF
	815
	816	! extract only the partitions which reduce the subdomain size in comparison with smaller partitions
	817	ALLOCATE( indexok(isz1) ) ! to store indices of the best partitions
[14072]	818	isz0 = 0 ! number of best partitions
[10425]	819	inbij = 1 ! start with the min value of inbij1 => 1
[13490]	820	iszij = jpiglo*jpjglo+1 ! default: larger than global domain
[10425]	821	DO WHILE( inbij <= inbijmax ) ! if we did not reach the max of inbij1
	822	ii = MINLOC(iszij1, mask = inbij1 == inbij, dim = 1) ! warning: send back the first occurence if multiple results
	823	IF ( iszij1(ii) < iszij ) THEN
[13490]	824	ii = MINLOC( iszi1+iszj1, mask = iszij1 == iszij1(ii) .AND. inbij1 == inbij, dim = 1) ! select the smaller perimeter if multiple min
[10425]	825	isz0 = isz0 + 1
	826	indexok(isz0) = ii
	827	iszij = iszij1(ii)
	828	ENDIF
	829	inbij = MINVAL(inbij1, mask = inbij1 > inbij) ! warning: return largest integer value if mask = .false. everywhere
	830	END DO
	831	DEALLOCATE( inbij1, iszij1 )
	832
	833	! keep only the best partitions (sorted by increasing order of subdomains number and decreassing subdomain size)
	834	ALLOCATE( inbi0(isz0), inbj0(isz0), iszi0(isz0), iszj0(isz0) )
	835	DO ji = 1, isz0
	836	ii = indexok(ji)
	837	inbi0(ji) = inbi1(ii)
	838	inbj0(ji) = inbj1(ii)
	839	iszi0(ji) = iszi1(ii)
	840	iszj0(ji) = iszj1(ii)
	841	END DO
	842	DEALLOCATE( indexok, inbi1, inbj1, iszi1, iszj1 )
	843
[11536]	844	IF( llist ) THEN
[10425]	845	IF(lwp) THEN
	846	WRITE(numout,*)
[11536]	847	WRITE(numout,*) ' For your information:'
	848	WRITE(numout,*) ' list of the best partitions including land supression'
	849	WRITE(numout,*) ' -----------------------------------------------------'
[10425]	850	WRITE(numout,*)
[14433]	851	ENDIF
[11536]	852	ji = isz0 ! initialization with the largest value
[14433]	853	ALLOCATE( llisOce(inbi0(ji), inbj0(ji)) )
	854	CALL mpp_is_ocean( llisOce ) ! Warning: must be call by all cores (call mpp_sum)
	855	inbijold = COUNT(llisOce)
	856	DEALLOCATE( llisOce )
[11536]	857	DO ji =isz0-1,1,-1
[14433]	858	ALLOCATE( llisOce(inbi0(ji), inbj0(ji)) )
	859	CALL mpp_is_ocean( llisOce ) ! Warning: must be call by all cores (call mpp_sum)
	860	inbij = COUNT(llisOce)
	861	DEALLOCATE( llisOce )
[11536]	862	IF(lwp .AND. inbij < inbijold) THEN
	863	WRITE(numout,'(a, i6, a, i6, a, f4.1, a, i9, a, i6, a, i6, a)') &
	864	& 'nb_cores oce: ', inbij, ', land domains excluded: ', inbi0(ji)*inbj0(ji) - inbij, &
	865	& ' (', REAL(inbi0(ji)inbj0(ji) - inbij,wp) / REAL(inbi0(ji)inbj0(ji),wp) *100., &
	866	& '%), largest oce domain: ', iszi0(ji)*iszj0(ji), ' ( ', iszi0(ji),' x ', iszj0(ji), ' )'
	867	inbijold = inbij
[14433]	868	ENDIF
[10425]	869	END DO
	870	DEALLOCATE( inbi0, inbj0, iszi0, iszj0 )
[11536]	871	IF(lwp) THEN
	872	WRITE(numout,*)
	873	WRITE(numout,*) ' -----------------------------------------------------------'
	874	ENDIF
	875	CALL mppsync
	876	CALL mppstop( ld_abort = .TRUE. )
[10425]	877	ENDIF
[14072]	878
[10425]	879	DEALLOCATE( iszi0, iszj0 )
	880	inbij = inbijmax + 1 ! default: larger than possible
	881	ii = isz0+1 ! start from the end of the list (smaller subdomains)
	882	DO WHILE( inbij > knbij ) ! while the number of ocean subdomains exceed the number of procs
[14072]	883	ii = ii -1
[14433]	884	ALLOCATE( llisOce(inbi0(ii), inbj0(ii)) )
	885	CALL mpp_is_ocean( llisOce ) ! must be done by all core
	886	inbij = COUNT(llisOce)
	887	DEALLOCATE( llisOce )
[10425]	888	END DO
	889	knbi = inbi0(ii)
	890	knbj = inbj0(ii)
	891	IF(PRESENT(knbcnt)) knbcnt = knbi * knbj - inbij
	892	DEALLOCATE( inbi0, inbj0 )
	893	!
[13286]	894	END SUBROUTINE bestpartition
[14072]	895
	896
[10425]	897	SUBROUTINE mpp_init_landprop( propland )
	898	!!----------------------------------------------------------------------
	899	!! * ROUTINE mpp_init_landprop *
	900	!!
	901	!! ** Purpose : the the proportion of land points in the surface land-sea mask
	902	!!
[13286]	903	!! ** Method : read iproc strips (of length Ni0glo) of the land-sea mask
[10425]	904	!!----------------------------------------------------------------------
[10436]	905	REAL(wp), INTENT( out) :: propland ! proportion of land points in the global domain (between 0 and 1)
[10425]	906	!
	907	INTEGER, DIMENSION(jpni*jpnj) :: kusedom_1d
[10436]	908	INTEGER :: inboce, iarea
[10425]	909	INTEGER :: iproc, idiv, ijsz
[10436]	910	INTEGER :: ijstr
[10425]	911	LOGICAL, ALLOCATABLE, DIMENSION(:,:) :: lloce
	912	!!----------------------------------------------------------------------
	913	! do nothing if there is no land-sea mask
	914	IF( numbot == -1 .and. numbdy == -1 ) THEN
	915	propland = 0.
	916	RETURN
	917	ENDIF
	918
[14072]	919	! number of processes reading the bathymetry file
[13286]	920	iproc = MINVAL( (/mppsize, Nj0glo/2, 100/) ) ! read a least 2 lines, no more that 100 processes reading at the same time
[14072]	921
[10436]	922	! we want to read iproc strips of the land-sea mask. -> pick up iproc processes every idiv processes starting at 1
[10425]	923	IF( iproc == 1 ) THEN ; idiv = mppsize
	924	ELSE ; idiv = ( mppsize - 1 ) / ( iproc - 1 )
	925	ENDIF
[10436]	926
	927	iarea = (narea-1)/idiv ! involed process number (starting counting at 0)
	928	IF( MOD( narea-1, idiv ) == 0 .AND. iarea < iproc ) THEN ! beware idiv can be = to 1
[10425]	929	!
[13286]	930	ijsz = Nj0glo / iproc ! width of the stripe to read
	931	IF( iarea < MOD(Nj0glo,iproc) ) ijsz = ijsz + 1
	932	ijstr = iarea*(Nj0glo/iproc) + MIN(iarea, MOD(Nj0glo,iproc)) + 1 ! starting j position of the reading
[10425]	933	!
[13286]	934	ALLOCATE( lloce(Ni0glo, ijsz) ) ! allocate the strip
[14433]	935	CALL read_mask( 1, ijstr, Ni0glo, ijsz, lloce )
[10436]	936	inboce = COUNT(lloce) ! number of ocean point in the stripe
[10425]	937	DEALLOCATE(lloce)
	938	!
	939	ELSE
	940	inboce = 0
	941	ENDIF
[10436]	942	CALL mpp_sum( 'mppini', inboce ) ! total number of ocean points over the global domain
[10425]	943	!
[14072]	944	propland = REAL( Ni0gloNj0glo - inboce, wp ) / REAL( Ni0gloNj0glo, wp )
[10425]	945	!
	946	END SUBROUTINE mpp_init_landprop
[14072]	947
	948
[14433]	949	SUBROUTINE mpp_is_ocean( ldIsOce )
[10425]	950	!!----------------------------------------------------------------------
[13286]	951	!! * ROUTINE mpp_is_ocean *
[10425]	952	!!
[13286]	953	!! ** Purpose : Check for a mpi domain decomposition inbi x inbj which
	954	!! subdomains, including 1 halo (even if nn_hls>1), contain
	955	!! at least 1 ocean point.
	956	!! We must indeed ensure that each subdomain that is a neighbour
[14433]	957	!! of a land subdomain, has only land points on its boundary
[13286]	958	!! (inside the inner subdomain) with the land subdomain.
	959	!! This is needed to get the proper bondary conditions on
	960	!! a subdomain with a closed boundary.
[10425]	961	!!
[13286]	962	!! ** Method : read inbj strips (of length Ni0glo) of the land-sea mask
[10425]	963	!!----------------------------------------------------------------------
[14433]	964	LOGICAL, DIMENSION(:,:), INTENT( out) :: ldIsOce ! .true. if a sub domain constains 1 ocean point
[10425]	965	!
[10436]	966	INTEGER :: idiv, iimax, ijmax, iarea
[13286]	967	INTEGER :: inbi, inbj, inx, iny, inry, isty
[10560]	968	INTEGER :: ji, jn
[13286]	969	INTEGER, ALLOCATABLE, DIMENSION(:,:) :: inboce ! number oce oce pint in each mpi subdomain
	970	INTEGER, ALLOCATABLE, DIMENSION(: ) :: inboce_1d
	971	INTEGER, ALLOCATABLE, DIMENSION(:,:) :: iimppt, ijpi
	972	INTEGER, ALLOCATABLE, DIMENSION(:,:) :: ijmppt, ijpj
[14072]	973	LOGICAL, ALLOCATABLE, DIMENSION(:,:) :: lloce ! lloce(i,j) = .true. if the point (i,j) is ocean
[10425]	974	!!----------------------------------------------------------------------
	975	! do nothing if there is no land-sea mask
	976	IF( numbot == -1 .AND. numbdy == -1 ) THEN
[14433]	977	ldIsOce(:,:) = .TRUE.
[10425]	978	RETURN
	979	ENDIF
[13286]	980	!
[14433]	981	inbi = SIZE( ldIsOce, dim = 1 )
	982	inbj = SIZE( ldIsOce, dim = 2 )
[13286]	983	!
	984	! we want to read inbj strips of the land-sea mask. -> pick up inbj processes every idiv processes starting at 1
	985	IF ( inbj == 1 ) THEN ; idiv = mppsize
	986	ELSE IF ( mppsize < inbj ) THEN ; idiv = 1
	987	ELSE ; idiv = ( mppsize - 1 ) / ( inbj - 1 )
[10425]	988	ENDIF
[13286]	989	!
	990	ALLOCATE( inboce(inbi,inbj), inboce_1d(inbi*inbj) )
[10436]	991	inboce(:,:) = 0 ! default no ocean point found
[13286]	992	!
	993	DO jn = 0, (inbj-1)/mppsize ! if mppsize < inbj : more strips than mpi processes (because of potential land domains)
[10425]	994	!
[13286]	995	iarea = (narea-1)/idiv + jn * mppsize + 1 ! involed process number (starting counting at 1)
	996	IF( MOD( narea-1, idiv ) == 0 .AND. iarea <= inbj ) THEN ! beware idiv can be = to 1
[10560]	997	!
[13286]	998	ALLOCATE( iimppt(inbi,inbj), ijmppt(inbi,inbj), ijpi(inbi,inbj), ijpj(inbi,inbj) )
	999	CALL mpp_basesplit( Ni0glo, Nj0glo, 0, inbi, inbj, iimax, ijmax, iimppt, ijmppt, ijpi, ijpj )
[10560]	1000	!
[13286]	1001	inx = Ni0glo + 2 ; iny = ijpj(1,iarea) + 2 ! strip size + 1 halo on each direction (even if nn_hls>1)
	1002	ALLOCATE( lloce(inx, iny) ) ! allocate the strip
	1003	inry = iny - COUNT( (/ iarea == 1, iarea == inbj /) ) ! number of point to read in y-direction
	1004	isty = 1 + COUNT( (/ iarea == 1 /) ) ! read from the first or the second line?
[14433]	1005	CALL read_mask( 1, ijmppt(1,iarea) - 2 + isty, Ni0glo, inry, lloce(2:inx-1, isty:inry+isty-1) ) ! read the strip
[14072]	1006	!
[13286]	1007	IF( iarea == 1 ) THEN ! the first line was not read
[14433]	1008	IF( l_Jperio ) THEN ! north-south periodocity
	1009	CALL read_mask( 1, Nj0glo, Ni0glo, 1, lloce(2:inx-1, 1) ) ! read the last line -> first line of lloce
[13286]	1010	ELSE
	1011	lloce(2:inx-1, 1) = .FALSE. ! closed boundary
	1012	ENDIF
	1013	ENDIF
	1014	IF( iarea == inbj ) THEN ! the last line was not read
[14433]	1015	IF( l_Jperio ) THEN ! north-south periodocity
	1016	CALL read_mask( 1, 1, Ni0glo, 1, lloce(2:inx-1,iny) ) ! read the first line -> last line of lloce
	1017	ELSEIF( c_NFtype == 'T' ) THEN ! north-pole folding T-pivot, T-point
[13286]	1018	lloce(2,iny) = lloce(2,iny-2) ! here we have 1 halo (even if nn_hls>1)
	1019	DO ji = 3,inx-1
	1020	lloce(ji,iny ) = lloce(inx-ji+2,iny-2) ! ok, we have at least 3 lines
	1021	END DO
	1022	DO ji = inx/2+2,inx-1
	1023	lloce(ji,iny-1) = lloce(inx-ji+2,iny-1)
	1024	END DO
[14433]	1025	ELSEIF( c_NFtype == 'F' ) THEN ! north-pole folding F-pivot, T-point, 1 halo
[13286]	1026	lloce(inx/2+1,iny-1) = lloce(inx/2,iny-1) ! here we have 1 halo (even if nn_hls>1)
	1027	lloce(inx -1,iny-1) = lloce(2 ,iny-1)
	1028	DO ji = 2,inx-1
	1029	lloce(ji,iny) = lloce(inx-ji+1,iny-1)
	1030	END DO
	1031	ELSE ! closed boundary
	1032	lloce(2:inx-1,iny) = .FALSE.
	1033	ENDIF
	1034	ENDIF
	1035	! ! first and last column were not read
[14433]	1036	IF( l_Iperio ) THEN
[13286]	1037	lloce(1,:) = lloce(inx-1,:) ; lloce(inx,:) = lloce(2,:) ! east-west periodocity
	1038	ELSE
	1039	lloce(1,:) = .FALSE. ; lloce(inx,:) = .FALSE. ! closed boundary
	1040	ENDIF
	1041	!
	1042	DO ji = 1, inbi
	1043	inboce(ji,iarea) = COUNT( lloce(iimppt(ji,1):iimppt(ji,1)+ijpi(ji,1)+1,:) ) ! lloce as 2 points more than Ni0glo
[10560]	1044	END DO
	1045	!
	1046	DEALLOCATE(lloce)
[13286]	1047	DEALLOCATE(iimppt, ijmppt, ijpi, ijpj)
[10560]	1048	!
	1049	ENDIF
	1050	END DO
[14072]	1051
[13286]	1052	inboce_1d = RESHAPE(inboce, (/ inbi*inbj /))
[10425]	1053	CALL mpp_sum( 'mppini', inboce_1d )
[13286]	1054	inboce = RESHAPE(inboce_1d, (/inbi, inbj/))
[14433]	1055	ldIsOce(:,:) = inboce(:,:) /= 0
[13286]	1056	DEALLOCATE(inboce, inboce_1d)
[10425]	1057	!
[13286]	1058	END SUBROUTINE mpp_is_ocean
[14072]	1059
	1060
[14433]	1061	SUBROUTINE read_mask( kistr, kjstr, kicnt, kjcnt, ldoce )
[10425]	1062	!!----------------------------------------------------------------------
[14433]	1063	!! * ROUTINE read_mask *
[10425]	1064	!!
	1065	!! ** Purpose : Read relevant bathymetric information in order to
	1066	!! provide a land/sea mask used for the elimination
	1067	!! of land domains, in an mpp computation.
	1068	!!
[13286]	1069	!! ** Method : read stipe of size (Ni0glo,...)
[10425]	1070	!!----------------------------------------------------------------------
[14433]	1071	INTEGER , INTENT(in ) :: kistr, kjstr ! starting i and j position of the reading
	1072	INTEGER , INTENT(in ) :: kicnt, kjcnt ! number of points to read in i and j directions
	1073	LOGICAL, DIMENSION(kicnt,kjcnt), INTENT( out) :: ldoce ! ldoce(i,j) = .true. if the point (i,j) is ocean
[10425]	1074	!
[14433]	1075	INTEGER :: inumsave ! local logical unit
	1076	REAL(wp), DIMENSION(kicnt,kjcnt) :: zbot, zbdy
[10425]	1077	!!----------------------------------------------------------------------
	1078	!
	1079	inumsave = numout ; numout = numnul ! redirect all print to /dev/null
	1080	!
[14072]	1081	IF( numbot /= -1 ) THEN
[14433]	1082	CALL iom_get( numbot, jpdom_unknown, 'bottom_level', zbot, kstart = (/kistr,kjstr/), kcount = (/kicnt, kjcnt/) )
[10425]	1083	ELSE
[13286]	1084	zbot(:,:) = 1._wp ! put a non-null value
[10425]	1085	ENDIF
[13286]	1086	!
[14072]	1087	IF( numbdy /= -1 ) THEN ! Adjust with bdy_msk if it exists
[14433]	1088	CALL iom_get ( numbdy, jpdom_unknown, 'bdy_msk', zbdy, kstart = (/kistr,kjstr/), kcount = (/kicnt, kjcnt/) )
[10425]	1089	zbot(:,:) = zbot(:,:) * zbdy(:,:)
	1090	ENDIF
	1091	!
[14433]	1092	ldoce(:,:) = NINT(zbot(:,:)) > 0
[10425]	1093	numout = inumsave
	1094	!
[14433]	1095	END SUBROUTINE read_mask
[10425]	1096
	1097
[14433]	1098	SUBROUTINE mpp_getnum( ldIsOce, kproc, kipos, kjpos )
[88]	1099	!!----------------------------------------------------------------------
[13286]	1100	!! * ROUTINE mpp_getnum *
[88]	1101	!!
[13286]	1102	!! ** Purpose : give a number to each MPI subdomains (starting at 0)
	1103	!!
	1104	!! ** Method : start from bottom left. First skip land subdomain, and finally use them if needed
	1105	!!----------------------------------------------------------------------
[14433]	1106	LOGICAL, DIMENSION(:,:), INTENT(in ) :: ldIsOce ! F if land process
	1107	INTEGER, DIMENSION(:,:), INTENT( out) :: kproc ! subdomain number (-1 if not existing, starting at 0)
[13286]	1108	INTEGER, DIMENSION( :), INTENT( out) :: kipos ! i-position of the subdomain (from 1 to jpni)
	1109	INTEGER, DIMENSION( :), INTENT( out) :: kjpos ! j-position of the subdomain (from 1 to jpnj)
	1110	!
	1111	INTEGER :: ii, ij, jarea, iarea0
	1112	INTEGER :: icont, i2add , ini, inj, inij
	1113	!!----------------------------------------------------------------------
	1114	!
[14433]	1115	ini = SIZE(ldIsOce, dim = 1)
	1116	inj = SIZE(ldIsOce, dim = 2)
[13286]	1117	inij = SIZE(kipos)
	1118	!
	1119	! specify which subdomains are oce subdomains; other are land subdomains
	1120	kproc(:,:) = -1
	1121	icont = -1
	1122	DO jarea = 1, ini*inj
	1123	iarea0 = jarea - 1
	1124	ii = 1 + MOD(iarea0,ini)
	1125	ij = 1 + iarea0/ini
[14433]	1126	IF( ldIsOce(ii,ij) ) THEN
[13286]	1127	icont = icont + 1
	1128	kproc(ii,ij) = icont
	1129	kipos(icont+1) = ii
	1130	kjpos(icont+1) = ij
	1131	ENDIF
	1132	END DO
	1133	! if needed add some land subdomains to reach inij active subdomains
[14433]	1134	i2add = inij - COUNT( ldIsOce )
[13286]	1135	DO jarea = 1, ini*inj
	1136	iarea0 = jarea - 1
	1137	ii = 1 + MOD(iarea0,ini)
	1138	ij = 1 + iarea0/ini
[14433]	1139	IF( .NOT. ldIsOce(ii,ij) .AND. i2add > 0 ) THEN
[13286]	1140	icont = icont + 1
	1141	kproc(ii,ij) = icont
	1142	kipos(icont+1) = ii
	1143	kjpos(icont+1) = ij
	1144	i2add = i2add - 1
	1145	ENDIF
	1146	END DO
	1147	!
	1148	END SUBROUTINE mpp_getnum
	1149
	1150
[14433]	1151	SUBROUTINE init_excl_landpt
	1152	!!----------------------------------------------------------------------
	1153	!! * ROUTINE *
	1154	!!
	1155	!! ** Purpose : exclude exchanges which contain only land points
	1156	!!
	1157	!! ** Method : if a send or receive buffer constains only land point we
	1158	!! flag off the corresponding communication
	1159	!! Warning: this selection depend on the halo size -> loop on halo size
	1160	!!
	1161	!!----------------------------------------------------------------------
	1162	INTEGER :: inumsave
	1163	INTEGER :: jh
	1164	INTEGER :: ipi, ipj
	1165	INTEGER :: iiwe, iiea, iist, iisz
	1166	INTEGER :: ijso, ijno, ijst, ijsz
	1167	LOGICAL :: llsave
	1168	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zmsk
	1169	LOGICAL , DIMENSION(Ni_0,Nj_0,1) :: lloce
	1170	!!----------------------------------------------------------------------
	1171	!
	1172	! read the land-sea mask on the inner domain
	1173	CALL read_mask( nimpp, njmpp, Ni_0, Nj_0, lloce(:,:,1) )
	1174	!
	1175	! Here we look only at communications excluding the NP folding.
	1176	! As lbcnfd not validated if halo size /= nn_hls, we switch if off temporary...
	1177	llsave = l_IdoNFold
	1178	l_IdoNFold = .FALSE.
	1179	!
	1180	DO jh = 1, n_hlsmax ! different halo size
	1181	!
	1182	ipi = Ni_0 + 2*jh ! local domain size
	1183	ipj = Nj_0 + 2*jh
	1184	!
	1185	ALLOCATE( zmsk(ipi,ipj) )
	1186	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
	1187	CALL lbc_lnk('mppini', zmsk, 'T', 1._wp, khls = jh) ! fill halos
	1188	!
	1189	iiwe = jh ; iiea = Ni_0 ! bottom-left corfer - 1 of the sent data
	1190	ijso = jh ; ijno = Nj_0
	1191	IF( nn_comm == 1 ) THEN
	1192	iist = 0 ; iisz = ipi
	1193	ijst = 0 ; ijsz = ipj
	1194	ELSE
	1195	iist = jh ; iisz = Ni_0
	1196	ijst = jh ; ijsz = Nj_0
	1197	ENDIF
	1198	IF( nn_comm == 1 ) THEN ! SM: NOT WORKING FOR NEIGHBOURHOOD COLLECTIVE COMMUNICATIONS, I DON'T KNOW WHY...
	1199	! do not send if we send only land points
	1200	IF( NINT(SUM( zmsk(iiwe+1:iiwe+jh ,ijst+1:ijst+ijsz) )) == 0 ) mpiSnei(jh,jpwe) = -1
	1201	IF( NINT(SUM( zmsk(iiea+1:iiea+jh ,ijst+1:ijst+ijsz) )) == 0 ) mpiSnei(jh,jpea) = -1
	1202	IF( NINT(SUM( zmsk(iist+1:iist+iisz,ijso+1:ijso+jh ) )) == 0 ) mpiSnei(jh,jpso) = -1
	1203	IF( NINT(SUM( zmsk(iist+1:iist+iisz,ijno+1:ijno+jh ) )) == 0 ) mpiSnei(jh,jpno) = -1
	1204	IF( NINT(SUM( zmsk(iiwe+1:iiwe+jh ,ijso+1:ijso+jh ) )) == 0 ) mpiSnei(jh,jpsw) = -1
	1205	IF( NINT(SUM( zmsk(iiea+1:iiea+jh ,ijso+1:ijso+jh ) )) == 0 ) mpiSnei(jh,jpse) = -1
	1206	IF( NINT(SUM( zmsk(iiwe+1:iiwe+jh ,ijno+1:ijno+jh ) )) == 0 ) mpiSnei(jh,jpnw) = -1
	1207	IF( NINT(SUM( zmsk(iiea+1:iiea+jh ,ijno+1:ijno+jh ) )) == 0 ) mpiSnei(jh,jpne) = -1
	1208	!
	1209	iiwe = iiwe-jh ; iiea = iiea+jh ! bottom-left corfer - 1 of the received data
	1210	ijso = ijso-jh ; ijno = ijno+jh
	1211	! do not send if we send only land points
	1212	IF( NINT(SUM( zmsk(iiwe+1:iiwe+jh ,ijst+1:ijst+ijsz) )) == 0 ) mpiRnei(jh,jpwe) = -1
	1213	IF( NINT(SUM( zmsk(iiea+1:iiea+jh ,ijst+1:ijst+ijsz) )) == 0 ) mpiRnei(jh,jpea) = -1
	1214	IF( NINT(SUM( zmsk(iist+1:iist+iisz,ijso+1:ijso+jh ) )) == 0 ) mpiRnei(jh,jpso) = -1
	1215	IF( NINT(SUM( zmsk(iist+1:iist+iisz,ijno+1:ijno+jh ) )) == 0 ) mpiRnei(jh,jpno) = -1
	1216	IF( NINT(SUM( zmsk(iiwe+1:iiwe+jh ,ijso+1:ijso+jh ) )) == 0 ) mpiRnei(jh,jpsw) = -1
	1217	IF( NINT(SUM( zmsk(iiea+1:iiea+jh ,ijso+1:ijso+jh ) )) == 0 ) mpiRnei(jh,jpse) = -1
	1218	IF( NINT(SUM( zmsk(iiwe+1:iiwe+jh ,ijno+1:ijno+jh ) )) == 0 ) mpiRnei(jh,jpnw) = -1
	1219	IF( NINT(SUM( zmsk(iiea+1:iiea+jh ,ijno+1:ijno+jh ) )) == 0 ) mpiRnei(jh,jpne) = -1
	1220	ENDIF
	1221	!
	1222	! Specific (and rare) problem in corner treatment because we do 1st West-East comm, next South-North comm
	1223	IF( nn_comm == 1 ) THEN
	1224	IF( mpiSnei(jh,jpwe) > -1 ) mpiSnei(jh, (/jpsw,jpnw/) ) = -1 ! SW and NW corners already sent through West nei
	1225	IF( mpiSnei(jh,jpea) > -1 ) mpiSnei(jh, (/jpse,jpne/) ) = -1 ! SE and NE corners already sent through East nei
	1226	IF( mpiRnei(jh,jpso) > -1 ) mpiRnei(jh, (/jpsw,jpse/) ) = -1 ! SW and SE corners will be received through South nei
	1227	IF( mpiRnei(jh,jpno) > -1 ) mpiRnei(jh, (/jpnw,jpne/) ) = -1 ! NW and NE corners will be received through North nei
	1228	ENDIF
	1229	!
	1230	DEALLOCATE( zmsk )
	1231	!
	1232	CALL mpp_ini_nc(jh) ! Initialize/Update communicator for neighbourhood collective communications
	1233	!
	1234	END DO
	1235	l_IdoNFold = llsave
	1236
	1237	END SUBROUTINE init_excl_landpt
	1238
	1239
[13286]	1240	SUBROUTINE init_ioipsl
	1241	!!----------------------------------------------------------------------
	1242	!! * ROUTINE init_ioipsl *
	1243	!!
[14072]	1244	!! ** Purpose :
[88]	1245	!!
[14072]	1246	!! ** Method :
[88]	1247	!!
	1248	!! History :
[14072]	1249	!! 9.0 ! 04-03 (G. Madec ) MPP-IOIPSL
[1238]	1250	!! " " ! 08-12 (A. Coward) addition in case of jpni*jpnj < jpnij
[88]	1251	!!----------------------------------------------------------------------
[2715]	1252	INTEGER, DIMENSION(2) :: iglo, iloc, iabsf, iabsl, ihals, ihale, idid
[88]	1253	!!----------------------------------------------------------------------
[352]	1254
[1238]	1255	! The domain is split only horizontally along i- or/and j- direction
	1256	! So we need at the most only 1D arrays with 2 elements.
	1257	! Set idompar values equivalent to the jpdom_local_noextra definition
	1258	! used in IOM. This works even if jpnij .ne. jpni*jpnj.
[13286]	1259	iglo( :) = (/ Ni0glo, Nj0glo /)
	1260	iloc( :) = (/ Ni_0 , Nj_0 /)
	1261	iabsf(:) = (/ Nis0 , Njs0 /) + (/ nimpp, njmpp /) - 1 - nn_hls ! corresponds to mig0(Nis0) but mig0 is not yet defined!
[88]	1262	iabsl(:) = iabsf(:) + iloc(:) - 1
[13286]	1263	ihals(:) = (/ 0 , 0 /)
	1264	ihale(:) = (/ 0 , 0 /)
	1265	idid( :) = (/ 1 , 2 /)
[352]	1266
[88]	1267	IF(lwp) THEN
[516]	1268	WRITE(numout,*)
[13286]	1269	WRITE(numout,*) 'mpp init_ioipsl : iloc = ', iloc
	1270	WRITE(numout,*) '~~~~~~~~~~~~~~~ iabsf = ', iabsf
	1271	WRITE(numout,*) ' ihals = ', ihals
	1272	WRITE(numout,*) ' ihale = ', ihale
[88]	1273	ENDIF
[2715]	1274	!
[14275]	1275	CALL flio_dom_set ( jpnij, narea-1, idid, iglo, iloc, iabsf, iabsl, ihals, ihale, 'BOX', nidom)
[2715]	1276	!
[14072]	1277	END SUBROUTINE init_ioipsl
[88]	1278
[9436]	1279
[13286]	1280	SUBROUTINE init_nfdcom
[9436]	1281	!!----------------------------------------------------------------------
[13286]	1282	!! * ROUTINE init_nfdcom *
[14072]	1283	!! ** Purpose : Setup for north fold exchanges with explicit
[9436]	1284	!! point-to-point messaging
	1285	!!
	1286	!! ** Method : Initialization of the northern neighbours lists.
	1287	!!----------------------------------------------------------------------
	1288	!! 1.0 ! 2011-10 (A. C. Coward, NOCS & J. Donners, PRACE)
[14072]	1289	!! 2.0 ! 2013-06 Setup avoiding MPI communication (I. Epicoco, S. Mocavero, CMCC)
[9436]	1290	!!----------------------------------------------------------------------
	1291	INTEGER :: sxM, dxM, sxT, dxT, jn
	1292	!!----------------------------------------------------------------------
	1293	!
[14433]	1294	!sxM is the first point (in the global domain) needed to compute the north-fold for the current process
	1295	sxM = jpiglo - nimpp - jpi + 1
	1296	!dxM is the last point (in the global domain) needed to compute the north-fold for the current process
	1297	dxM = jpiglo - nimpp + 2
[9436]	1298	!
[14433]	1299	! loop over the other north-fold processes to find the processes
	1300	! managing the points belonging to the sxT-dxT range
	1301	!
	1302	nsndto = 0
	1303	DO jn = 1, jpni
[9436]	1304	!
[14433]	1305	sxT = nfimpp(jn) ! sxT = 1st point (in the global domain) of the jn process
	1306	dxT = nfimpp(jn) + nfjpi(jn) - 1 ! dxT = last point (in the global domain) of the jn process
[9436]	1307	!
[14433]	1308	IF ( sxT < sxM .AND. sxM < dxT ) THEN
	1309	nsndto = nsndto + 1
	1310	isendto(nsndto) = jn
	1311	ELSEIF( sxM <= sxT .AND. dxM >= dxT ) THEN
	1312	nsndto = nsndto + 1
	1313	isendto(nsndto) = jn
	1314	ELSEIF( dxM < dxT .AND. sxT < dxM ) THEN
	1315	nsndto = nsndto + 1
	1316	isendto(nsndto) = jn
	1317	ENDIF
[9436]	1318	!
[14433]	1319	END DO
[9436]	1320	!
[13286]	1321	END SUBROUTINE init_nfdcom
[9436]	1322
[88]	1323
[13286]	1324	SUBROUTINE init_doloop
	1325	!!----------------------------------------------------------------------
	1326	!! * ROUTINE init_doloop *
	1327	!!
	1328	!! ** Purpose : set the starting/ending indices of DO-loop
	1329	!! These indices are used in do_loop_substitute.h90
	1330	!!----------------------------------------------------------------------
	1331	!
[14433]	1332	Nis0 = 1+nn_hls
	1333	Njs0 = 1+nn_hls
	1334	Nie0 = jpi-nn_hls
	1335	Nje0 = jpj-nn_hls
[14072]	1336	!
[13286]	1337	Ni_0 = Nie0 - Nis0 + 1
	1338	Nj_0 = Nje0 - Njs0 + 1
	1339	!
[14433]	1340	jpkm1 = jpk-1 ! " "
	1341	!
[13286]	1342	END SUBROUTINE init_doloop
[14072]	1343
[3]	1344	!!======================================================================
	1345	END MODULE mppini

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

Download in other formats: