MODULE mppini !!====================================================================== !! *** MODULE mppini *** !! Ocean initialization : distributed memory computing initialization !!====================================================================== !! History : 6.0 ! 1994-11 (M. Guyon) Original code !! OPA 7.0 ! 1995-04 (J. Escobar, M. Imbard) !! 8.0 ! 1998-05 (M. Imbard, J. Escobar, L. Colombet ) SHMEM and MPI versions !! NEMO 1.0 ! 2004-01 (G. Madec, J.M Molines) F90 : free form , north fold jpni > 1 !! 3.4 ! 2011-10 (A. C. Coward, NOCS & J. Donners, PRACE) add init_nfdcom !! 3. ! 2013-06 (I. Epicoco, S. Mocavero, CMCC) init_nfdcom: setup avoiding MPI communication !! 4.0 ! 2016-06 (G. Madec) use domain configuration file instead of bathymetry file !! 4.0 ! 2017-06 (J.M. Molines, T. Lovato) merge of mppini and mppini_2 !!---------------------------------------------------------------------- !!---------------------------------------------------------------------- !! mpp_init : Lay out the global domain over processors with/without land processor elimination !! init_ioipsl: IOIPSL initialization in mpp !! init_nfdcom: Setup for north fold exchanges with explicit point-to-point messaging !! init_doloop: set the starting/ending indices of DO-loop used in do_loop_substitute !!---------------------------------------------------------------------- USE dom_oce ! ocean space and time domain USE bdy_oce ! open BounDarY ! USE lbcnfd , ONLY : isendto, nsndto ! Setup of north fold exchanges USE lib_mpp ! distribued memory computing library USE iom ! nemo I/O library USE ioipsl ! I/O IPSL library USE in_out_manager ! I/O Manager IMPLICIT NONE PRIVATE PUBLIC mpp_init ! called by nemogcm.F90 PUBLIC mpp_getnum ! called by prtctl PUBLIC mpp_basesplit ! called by prtctl PUBLIC mpp_is_ocean ! called by prtctl INTEGER :: numbot = -1 ! 'bottom_level' local logical unit INTEGER :: numbdy = -1 ! 'bdy_msk' local logical unit !!---------------------------------------------------------------------- !! NEMO/OCE 4.0 , NEMO Consortium (2018) !! $Id$ !! Software governed by the CeCILL license (see ./LICENSE) !!---------------------------------------------------------------------- CONTAINS #if defined key_mpi_off !!---------------------------------------------------------------------- !! Default option : shared memory computing !!---------------------------------------------------------------------- SUBROUTINE mpp_init !!---------------------------------------------------------------------- !! *** ROUTINE mpp_init *** !! !! ** Purpose : Lay out the global domain over processors. !! !! ** Method : Shared memory computing, set the local processor !! variables to the value of the global domain !!---------------------------------------------------------------------- ! nn_hls = 1 jpiglo = Ni0glo + 2 * nn_hls jpjglo = Nj0glo + 2 * nn_hls jpimax = jpiglo jpjmax = jpjglo jpi = jpiglo jpj = jpjglo jpk = MAX( 2, jpkglo ) jpij = jpi*jpj jpni = 1 jpnj = 1 jpnij = jpni*jpnj nimpp = 1 njmpp = 1 nidom = FLIO_DOM_NONE ! CALL init_doloop ! set start/end indices or do-loop depending on the halo width value (nn_hls) ! IF(lwp) THEN WRITE(numout,*) WRITE(numout,*) 'mpp_init : NO massively parallel processing' WRITE(numout,*) '~~~~~~~~ ' WRITE(numout,*) ' l_Iperio = ', l_Iperio, ' l_Jperio = ', l_Jperio WRITE(numout,*) ' njmpp = ', njmpp ENDIF ! #if defined key_agrif IF (.NOT.agrif_root()) THEN call agrif_nemo_init() ENDIF #endif END SUBROUTINE mpp_init #else !!---------------------------------------------------------------------- !! MPI massively parallel processing !!---------------------------------------------------------------------- SUBROUTINE mpp_init !!---------------------------------------------------------------------- !! *** ROUTINE mpp_init *** !! !! ** Purpose : Lay out the global domain over processors. !! If land processors are to be eliminated, this program requires the !! presence of the domain configuration file. Land processors elimination !! is performed if jpni x jpnj /= jpnij. In this case, using the MPP_PREP !! preprocessing tool, help for defining the best cutting out. !! !! ** Method : Global domain is distributed in smaller local domains. !! Periodic condition is a function of the local domain position !! (global boundary or neighbouring domain) and of the global periodic !! !! ** Action : - set domain parameters !! nimpp : longitudinal index !! njmpp : latitudinal index !! narea : number for local area !! mpinei : number of neighboring domains (starting at 0, -1 if no neighbourg) !!---------------------------------------------------------------------- INTEGER :: ji, jj, jn, jp, jh INTEGER :: ii, ij, ii2, ij2 INTEGER :: inijmin ! number of oce subdomains INTEGER :: inum, inum0 INTEGER :: ifreq, il1, imil, il2, ijm1 INTEGER :: ierr, ios INTEGER :: inbi, inbj, iimax, ijmax, icnt1, icnt2 INTEGER, DIMENSION(16*n_hlsmax) :: ichanged INTEGER, ALLOCATABLE, DIMENSION(: ) :: iin, ijn INTEGER, ALLOCATABLE, DIMENSION(:,: ) :: iimppt, ijpi, ipproc INTEGER, ALLOCATABLE, DIMENSION(:,: ) :: ijmppt, ijpj INTEGER, ALLOCATABLE, DIMENSION(:,: ) :: impi INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: inei LOGICAL :: llbest, llauto LOGICAL :: llwrtlay LOGICAL :: llmpi_Iperio, llmpi_Jperio, llmpiNFold LOGICAL :: ln_listonly LOGICAL, ALLOCATABLE, DIMENSION(:,: ) :: llisOce ! is not land-domain only? LOGICAL, ALLOCATABLE, DIMENSION(:,:,:) :: llnei ! are neighbourgs existing? NAMELIST/nambdy/ ln_bdy, nb_bdy, ln_coords_file, cn_coords_file, & & ln_mask_file, cn_mask_file, cn_dyn2d, nn_dyn2d_dta, & & cn_dyn3d, nn_dyn3d_dta, cn_tra, nn_tra_dta, & & ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp, rn_time_dmp_out, & & cn_ice, nn_ice_dta, & & ln_vol, nn_volctl, nn_rimwidth NAMELIST/nammpp/ jpni, jpnj, nn_hls, ln_nnogather, ln_listonly, nn_comm !!---------------------------------------------------------------------- ! llwrtlay = lwm .OR. sn_cfctl%l_layout ! ! 0. read namelists parameters ! ----------------------------------- ! READ ( numnam_ref, nammpp, IOSTAT = ios, ERR = 901 ) 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nammpp in reference namelist' ) READ ( numnam_cfg, nammpp, IOSTAT = ios, ERR = 902 ) 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'nammpp in configuration namelist' ) ! nn_hls = MAX(1, nn_hls) ! nn_hls must be > 0 IF(lwp) THEN WRITE(numout,*) ' Namelist nammpp' IF( jpni < 1 .OR. jpnj < 1 ) THEN WRITE(numout,*) ' jpni and jpnj will be calculated automatically' ELSE WRITE(numout,*) ' processor grid extent in i jpni = ', jpni WRITE(numout,*) ' processor grid extent in j jpnj = ', jpnj ENDIF WRITE(numout,*) ' avoid use of mpi_allgather at the north fold ln_nnogather = ', ln_nnogather WRITE(numout,*) ' halo width (applies to both rows and columns) nn_hls = ', nn_hls ENDIF ! IF(lwm) WRITE( numond, nammpp ) ! jpiglo = Ni0glo + 2 * nn_hls jpjglo = Nj0glo + 2 * nn_hls ! ! do we need to take into account bdy_msk? READ ( numnam_ref, nambdy, IOSTAT = ios, ERR = 903) 903 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy in reference namelist (mppini)' ) READ ( numnam_cfg, nambdy, IOSTAT = ios, ERR = 904 ) 904 IF( ios > 0 ) CALL ctl_nam ( ios , 'nambdy in configuration namelist (mppini)' ) ! IF( ln_read_cfg ) CALL iom_open( cn_domcfg, numbot ) IF( ln_bdy .AND. ln_mask_file ) CALL iom_open( cn_mask_file, numbdy ) ! IF( ln_listonly ) CALL bestpartition( MAX(mppsize,jpni*jpnj), ldlist = .TRUE. ) ! must be done by all core ! ! 1. Dimension arrays for subdomains ! ----------------------------------- ! ! If dimensions of MPI processes grid weren't specified in the namelist file ! then we calculate them here now that we have our communicator size IF(lwp) THEN WRITE(numout,*) WRITE(numout,*) 'mpp_init:' WRITE(numout,*) '~~~~~~~~ ' ENDIF IF( jpni < 1 .OR. jpnj < 1 ) THEN CALL bestpartition( mppsize, jpni, jpnj ) ! best mpi decomposition for mppsize mpi processes llauto = .TRUE. llbest = .TRUE. ELSE llauto = .FALSE. CALL bestpartition( mppsize, inbi, inbj, icnt2 ) ! best mpi decomposition for mppsize mpi processes ! largest subdomain size for mpi decoposition jpni*jpnj given in the namelist CALL mpp_basesplit( jpiglo, jpjglo, nn_hls, jpni, jpnj, jpimax, jpjmax ) ! largest subdomain size for mpi decoposition inbi*inbj given by bestpartition CALL mpp_basesplit( jpiglo, jpjglo, nn_hls, inbi, inbj, iimax, ijmax ) icnt1 = jpni*jpnj - mppsize ! number of land subdomains that should be removed to use mppsize mpi processes IF(lwp) THEN WRITE(numout,9000) ' The chosen domain decomposition ', jpni, ' x ', jpnj, ' with ', icnt1, ' land subdomains' WRITE(numout,9002) ' - uses a total of ', mppsize,' mpi process' WRITE(numout,9000) ' - has mpi subdomains with a maximum size of (jpi = ', jpimax, ', jpj = ', jpjmax, & & ', jpi*jpj = ', jpimax*jpjmax, ')' WRITE(numout,9000) ' The best domain decompostion ', inbi, ' x ', inbj, ' with ', icnt2, ' land subdomains' WRITE(numout,9002) ' - uses a total of ', inbi*inbj-icnt2,' mpi process' WRITE(numout,9000) ' - has mpi subdomains with a maximum size of (jpi = ', iimax, ', jpj = ', ijmax, & & ', jpi*jpj = ', iimax* ijmax, ')' ENDIF IF( iimax*ijmax < jpimax*jpjmax ) THEN ! chosen subdomain size is larger that the best subdomain size llbest = .FALSE. IF ( inbi*inbj-icnt2 < mppsize ) THEN WRITE(ctmp1,*) ' ==> You could therefore have smaller mpi subdomains with less mpi processes' ELSE WRITE(ctmp1,*) ' ==> You could therefore have smaller mpi subdomains with the same number of mpi processes' ENDIF CALL ctl_warn( ' ', ctmp1, ' ', ' --- YOU ARE WASTING CPU... ---', ' ' ) ELSE IF ( iimax*ijmax == jpimax*jpjmax .AND. (inbi*inbj-icnt2) < mppsize) THEN llbest = .FALSE. WRITE(ctmp1,*) ' ==> You could therefore have the same mpi subdomains size with less mpi processes' CALL ctl_warn( ' ', ctmp1, ' ', ' --- YOU ARE WASTING CPU... ---', ' ' ) ELSE llbest = .TRUE. ENDIF ENDIF ! look for land mpi subdomains... ALLOCATE( llisOce(jpni,jpnj) ) CALL mpp_is_ocean( llisOce ) inijmin = COUNT( llisOce ) ! number of oce subdomains IF( mppsize < inijmin ) THEN ! too many oce subdomains: can happen only if jpni and jpnj are prescribed... WRITE(ctmp1,9001) ' With this specified domain decomposition: jpni = ', jpni, ' jpnj = ', jpnj WRITE(ctmp2,9002) ' we can eliminate only ', jpni*jpnj - inijmin, ' land mpi subdomains therefore ' WRITE(ctmp3,9001) ' the number of ocean mpi subdomains (', inijmin,') exceed the number of MPI processes:', mppsize WRITE(ctmp4,*) ' ==>>> There is the list of best domain decompositions you should use: ' CALL ctl_stop( ctmp1, ctmp2, ctmp3, ' ', ctmp4, ' ' ) CALL bestpartition( mppsize, ldlist = .TRUE. ) ! must be done by all core ENDIF IF( mppsize > jpni*jpnj ) THEN ! not enough mpi subdomains for the total number of mpi processes IF(lwp) THEN WRITE(numout,9003) ' The number of mpi processes: ', mppsize WRITE(numout,9003) ' exceeds the maximum number of subdomains (ocean+land) = ', jpni*jpnj WRITE(numout,9001) ' defined by the following domain decomposition: jpni = ', jpni, ' jpnj = ', jpnj WRITE(numout, *) ' You should: ' IF( llauto ) THEN WRITE(numout,*) ' - either prescribe your domain decomposition with the namelist variables' WRITE(numout,*) ' jpni and jpnj to match the number of mpi process you want to use, ' WRITE(numout,*) ' even IF it not the best choice...' WRITE(numout,*) ' - or keep the automatic and optimal domain decomposition by picking up one' WRITE(numout,*) ' of the number of mpi process proposed in the list bellow' ELSE WRITE(numout,*) ' - either properly prescribe your domain decomposition with jpni and jpnj' WRITE(numout,*) ' in order to be consistent with the number of mpi process you want to use' WRITE(numout,*) ' even IF it not the best choice...' WRITE(numout,*) ' - or use the automatic and optimal domain decomposition and pick up one of' WRITE(numout,*) ' the domain decomposition proposed in the list bellow' ENDIF WRITE(numout,*) ENDIF CALL bestpartition( mppsize, ldlist = .TRUE. ) ! must be done by all core ENDIF jpnij = mppsize ! force jpnij definition <-- remove as much land subdomains as needed to reach this condition IF( mppsize > inijmin ) THEN WRITE(ctmp1,9003) ' The number of mpi processes: ', mppsize WRITE(ctmp2,9003) ' exceeds the maximum number of ocean subdomains = ', inijmin WRITE(ctmp3,9002) ' we suppressed ', jpni*jpnj - mppsize, ' land subdomains ' WRITE(ctmp4,9002) ' BUT we had to keep ', mppsize - inijmin, ' land subdomains that are useless...' CALL ctl_warn( ctmp1, ctmp2, ctmp3, ctmp4, ' ', ' --- YOU ARE WASTING CPU... ---', ' ' ) ELSE ! mppsize = inijmin IF(lwp) THEN IF(llbest) WRITE(numout,*) ' ==> you use the best mpi decomposition' WRITE(numout,*) WRITE(numout,9003) ' Number of mpi processes: ', mppsize WRITE(numout,9003) ' Number of ocean subdomains = ', inijmin WRITE(numout,9003) ' Number of suppressed land subdomains = ', jpni*jpnj - inijmin WRITE(numout,*) ENDIF ENDIF 9000 FORMAT (a, i4, a, i4, a, i7, a) 9001 FORMAT (a, i4, a, i4) 9002 FORMAT (a, i4, a) 9003 FORMAT (a, i5) ALLOCATE( nfimpp(jpni), nfproc(jpni), nfjpi(jpni), & & iin(jpnij), ijn(jpnij), & & iimppt(jpni,jpnj), ijmppt(jpni,jpnj), ijpi(jpni,jpnj), ijpj(jpni,jpnj), ipproc(jpni,jpnj), & & inei(8,jpni,jpnj), llnei(8,jpni,jpnj), & & impi(8,jpnij), & & STAT=ierr ) CALL mpp_sum( 'mppini', ierr ) IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'mpp_init: unable to allocate standard ocean arrays' ) #if defined key_agrif IF( .NOT. Agrif_Root() ) THEN ! AGRIF children: specific setting (cf. agrif_user.F90) CALL agrif_nemo_init() ENDIF #endif ! ! 2. Index arrays for subdomains ! ----------------------------------- ! CALL mpp_basesplit( jpiglo, jpjglo, nn_hls, jpni, jpnj, jpimax, jpjmax, iimppt, ijmppt, ijpi, ijpj ) CALL mpp_getnum( llisOce, ipproc, iin, ijn ) ! ii = iin(narea) ij = ijn(narea) jpi = ijpi(ii,ij) jpj = ijpj(ii,ij) jpk = MAX( 2, jpkglo ) jpij = jpi*jpj nimpp = iimppt(ii,ij) njmpp = ijmppt(ii,ij) ! CALL init_doloop ! set start/end indices of do-loop, depending on the halo width value (nn_hls) ! IF(lwp) THEN WRITE(numout,*) WRITE(numout,*) 'MPI Message Passing MPI - domain lay out over processors' WRITE(numout,*) WRITE(numout,*) ' defines mpp subdomains' WRITE(numout,*) ' jpni = ', jpni WRITE(numout,*) ' jpnj = ', jpnj WRITE(numout,*) ' jpnij = ', jpnij WRITE(numout,*) ' nimpp = ', nimpp WRITE(numout,*) ' njmpp = ', njmpp WRITE(numout,*) WRITE(numout,*) ' sum ijpi(i,1) = ', sum(ijpi(:,1)), ' jpiglo = ', jpiglo WRITE(numout,*) ' sum ijpj(1,j) = ', SUM(ijpj(1,:)), ' jpjglo = ', jpjglo ! Subdomain grid print ifreq = 4 il1 = 1 DO jn = 1, (jpni-1)/ifreq+1 il2 = MIN(jpni,il1+ifreq-1) WRITE(numout,*) WRITE(numout,9400) ('***',ji=il1,il2-1) DO jj = jpnj, 1, -1 WRITE(numout,9403) (' ',ji=il1,il2-1) WRITE(numout,9402) jj, (ijpi(ji,jj),ijpj(ji,jj),ji=il1,il2) WRITE(numout,9404) (ipproc(ji,jj),ji=il1,il2) WRITE(numout,9403) (' ',ji=il1,il2-1) WRITE(numout,9400) ('***',ji=il1,il2-1) END DO WRITE(numout,9401) (ji,ji=il1,il2) il1 = il1+ifreq END DO 9400 FORMAT(' ***' ,20('*************',a3) ) 9403 FORMAT(' * ',20(' * ',a3) ) 9401 FORMAT(' ' ,20(' ',i3,' ') ) 9402 FORMAT(' ',i3,' * ',20(i3,' x',i3,' * ') ) 9404 FORMAT(' * ' ,20(' ' ,i4,' * ') ) ENDIF ! ! Store informations for the north pole folding communications nfproc(:) = ipproc(:,jpnj) nfimpp(:) = iimppt(:,jpnj) nfjpi (:) = ijpi(:,jpnj) ! ! 3. Define Western, Eastern, Southern and Northern neighbors + corners in the subdomain grid reference ! ------------------------------------------------------------------------------------------------------ ! ! note that North fold is has specific treatment for its MPI communications. ! This must not be treated as a "usual" communication with a northern neighbor. ! -> North fold processes have no Northern neighbor in the definition done bellow ! llmpi_Iperio = jpni > 1 .AND. l_Iperio ! do i-periodicity with an MPI communication? llmpi_Jperio = jpnj > 1 .AND. l_Jperio ! do j-periodicity with an MPI communication? ! l_SelfPerio(1:2) = l_Iperio .AND. jpni == 1 ! west, east periodicity by itself l_SelfPerio(3:4) = l_Jperio .AND. jpnj == 1 ! south, north periodicity by itself l_SelfPerio(5:8) = l_SelfPerio(jpwe) .AND. l_SelfPerio(jpso) ! corners bi-periodicity by itself ! ! define neighbors mapping (1/2): default definition: ignore if neighbours are land-only subdomains or not DO jj = 1, jpnj DO ji = 1, jpni ! IF ( llisOce(ji,jj) ) THEN ! this subdomain has some ocean: it has neighbours ! inum0 = ji - 1 + ( jj - 1 ) * jpni ! index in the subdomains grid. start at 0 ! ! Is there a neighbor? llnei(jpwe,ji,jj) = ji > 1 .OR. llmpi_Iperio ! West nei exists if not the first column or llmpi_Iperio llnei(jpea,ji,jj) = ji < jpni .OR. llmpi_Iperio ! East nei exists if not the last column or llmpi_Iperio llnei(jpso,ji,jj) = jj > 1 .OR. llmpi_Jperio ! South nei exists if not the first line or llmpi_Jperio llnei(jpno,ji,jj) = jj < jpnj .OR. llmpi_Jperio ! North nei exists if not the last line or llmpi_Jperio llnei(jpsw,ji,jj) = llnei(jpwe,ji,jj) .AND. llnei(jpso,ji,jj) ! So-We nei exists if both South and West nei exist llnei(jpse,ji,jj) = llnei(jpea,ji,jj) .AND. llnei(jpso,ji,jj) ! So-Ea nei exists if both South and East nei exist llnei(jpnw,ji,jj) = llnei(jpwe,ji,jj) .AND. llnei(jpno,ji,jj) ! No-We nei exists if both North and West nei exist llnei(jpne,ji,jj) = llnei(jpea,ji,jj) .AND. llnei(jpno,ji,jj) ! No-Ea nei exists if both North and East nei exist ! ! Which index (starting at 0) have neighbors in the subdomains grid? IF( llnei(jpwe,ji,jj) ) inei(jpwe,ji,jj) = inum0 - 1 + jpni * COUNT( (/ ji == 1 /) ) IF( llnei(jpea,ji,jj) ) inei(jpea,ji,jj) = inum0 + 1 - jpni * COUNT( (/ ji == jpni /) ) IF( llnei(jpso,ji,jj) ) inei(jpso,ji,jj) = inum0 - jpni + jpni * jpnj * COUNT( (/ jj == 1 /) ) IF( llnei(jpno,ji,jj) ) inei(jpno,ji,jj) = inum0 + jpni - jpni * jpnj * COUNT( (/ jj == jpnj /) ) IF( llnei(jpsw,ji,jj) ) inei(jpsw,ji,jj) = inei(jpso,ji,jj) - 1 + jpni * COUNT( (/ ji == 1 /) ) IF( llnei(jpse,ji,jj) ) inei(jpse,ji,jj) = inei(jpso,ji,jj) + 1 - jpni * COUNT( (/ ji == jpni /) ) IF( llnei(jpnw,ji,jj) ) inei(jpnw,ji,jj) = inei(jpno,ji,jj) - 1 + jpni * COUNT( (/ ji == 1 /) ) IF( llnei(jpne,ji,jj) ) inei(jpne,ji,jj) = inei(jpno,ji,jj) + 1 - jpni * COUNT( (/ ji == jpni /) ) ! ELSE ! land-only domain has no neighbour llnei(:,ji,jj) = .FALSE. ENDIF ! END DO END DO ! ! define neighbors mapping (2/2): check if neighbours are not land-only subdomains DO jj = 1, jpnj DO ji = 1, jpni DO jn = 1, 8 IF( llnei(jn,ji,jj) ) THEN ! if a neighbour is existing -> this should not be a land-only domain ii = 1 + MOD( inei(jn,ji,jj) , jpni ) ij = 1 + inei(jn,ji,jj) / jpni llnei(jn,ji,jj) = llisOce( ii, ij ) ENDIF END DO END DO END DO ! ! update index of the neighbours in the subdomains grid WHERE( .NOT. llnei ) inei = -1 ! ! Save processor layout in ascii file IF (llwrtlay) THEN CALL ctl_opn( inum, 'layout.dat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE., narea ) WRITE(inum,'(a)') ' jpnij jpimax jpjmax jpk jpiglo jpjglo ( local: narea jpi jpj )' WRITE(inum,'(6i7,a,3i7,a)') jpnij,jpimax,jpjmax,jpk,jpiglo,jpjglo,' ( local: ',narea,jpi,jpj,' )' WRITE(inum,*) WRITE(inum, *) '------------------------------------' WRITE(inum,'(a,i2)') ' Mapping of the default neighnourgs ' WRITE(inum, *) '------------------------------------' WRITE(inum,*) WRITE(inum,'(a)') ' rank ii ij jpi jpj nimpp njmpp mpiwe mpiea mpiso mpino mpisw mpise mpinw mpine' DO jp = 1, jpnij ii = iin(jp) ij = ijn(jp) WRITE(inum,'(15i6)') jp-1, ii, ij, ijpi(ii,ij), ijpj(ii,ij), iimppt(ii,ij), ijmppt(ii,ij), inei(:,ii,ij) END DO ENDIF ! ! 4. Define Western, Eastern, Southern and Northern neighbors + corners for each mpi process ! ------------------------------------------------------------------------------------------ ! ! rewrite information from "subdomain grid" to mpi process list ! Warning, for example: ! position of the northern neighbor in the "subdomain grid" ! position of the northern neighbor in the "mpi process list" ! default definition: no neighbors impi(:,:) = -1 ! (starting at 0, -1 if no neighbourg) DO jp = 1, jpnij ii = iin(jp) ij = ijn(jp) DO jn = 1, 8 IF( llnei(jn,ii,ij) ) THEN ! must be tested as some land-domain can be kept to fit mppsize ii2 = 1 + MOD( inei(jn,ii,ij) , jpni ) ij2 = 1 + inei(jn,ii,ij) / jpni impi(jn,jp) = ipproc( ii2, ij2 ) ENDIF END DO END DO ! ! 4. keep information for the local process ! ----------------------------------------- ! ! set default neighbours mpinei(:) = impi(:,narea) DO jh = 1, n_hlsmax mpiSnei(jh,:) = impi(:,narea) ! default definition mpiRnei(jh,:) = impi(:,narea) END DO ! IF(lwp) THEN WRITE(numout,*) WRITE(numout,*) ' resulting internal parameters : ' WRITE(numout,*) ' narea = ', narea WRITE(numout,*) ' mpi nei west = ', mpinei(jpwe) , ' mpi nei east = ', mpinei(jpea) WRITE(numout,*) ' mpi nei south = ', mpinei(jpso) , ' mpi nei north = ', mpinei(jpno) WRITE(numout,*) ' mpi nei so-we = ', mpinei(jpsw) , ' mpi nei so-ea = ', mpinei(jpse) WRITE(numout,*) ' mpi nei no-we = ', mpinei(jpnw) , ' mpi nei no-ea = ', mpinei(jpne) ENDIF ! ! Prepare mpp north fold ! llmpiNFold = jpni > 1 .AND. l_NFold ! is the North fold done with an MPI communication? l_IdoNFold = ijn(narea) == jpnj .AND. l_NFold ! is this process doing North fold? ! IF( llmpiNFold ) THEN CALL mpp_ini_north IF (lwp) THEN WRITE(numout,*) WRITE(numout,*) ' ==>>> North fold boundary prepared for jpni >1' ENDIF IF (llwrtlay) THEN ! additional prints in layout.dat WRITE(inum,*) WRITE(inum,*) WRITE(inum,*) 'Number of subdomains located along the north fold : ', ndim_rank_north WRITE(inum,*) 'Rank of the subdomains located along the north fold : ', ndim_rank_north DO jp = 1, ndim_rank_north, 5 WRITE(inum,*) nrank_north( jp:MINVAL( (/jp+4,ndim_rank_north/) ) ) END DO ENDIF IF ( l_IdoNFold .AND. ln_nnogather ) THEN CALL init_nfdcom ! northfold neighbour lists IF (llwrtlay) THEN WRITE(inum,*) WRITE(inum,*) 'north fold exchanges with explicit point-to-point messaging :' WRITE(inum,*) ' nsndto : ', nsndto WRITE(inum,*) ' isendto : ', isendto(1:nsndto) ENDIF ENDIF ENDIF ! CALL mpp_ini_nc(nn_hls) ! Initialize communicator for neighbourhood collective communications DO jh = 1, n_hlsmax mpi_nc_com4(jh) = mpi_nc_com4(nn_hls) ! default definition mpi_nc_com8(jh) = mpi_nc_com8(nn_hls) END DO ! CALL init_excl_landpt ! exclude exchanges which contain only land points ! ! Save processor layout changes in ascii file DO jh = 1, n_hlsmax ! different halo size DO ji = 1, 8 ichanged(16*(jh-1) +ji) = COUNT( mpinei(ji:ji) /= mpiSnei(jh,ji:ji) ) ichanged(16*(jh-1)+8+ji) = COUNT( mpinei(ji:ji) /= mpiRnei(jh,ji:ji) ) END DO END DO CALL mpp_sum( "mpp_init", ichanged ) ! must be called by all processes IF (llwrtlay) THEN WRITE(inum,*) WRITE(inum, *) '----------------------------------------------------------------------' WRITE(inum,'(a,i2)') ' Mapping of the neighnourgs once excluding comm with only land points ' WRITE(inum, *) '----------------------------------------------------------------------' DO jh = 1, n_hlsmax ! different halo size WRITE(inum,*) WRITE(inum,'(a,i2)') 'halo size: ', jh WRITE(inum, *) '---------' WRITE(inum,'(a)') ' rank ii ij mpiwe mpiea mpiso mpino mpisw mpise mpinw mpine' WRITE(inum, '(11i6,a)') narea-1, iin(narea), ijn(narea), mpinei(:), ' <- Org' WRITE(inum,'(18x,8i6,a,i1,a)') mpiSnei(jh,:), ' <- Send ', COUNT( mpinei(:) /= mpiSnei(jh,:) ), ' modif' WRITE(inum,'(18x,8i6,a,i1,a)') mpiRnei(jh,:), ' <- Recv ', COUNT( mpinei(:) /= mpiRnei(jh,:) ), ' modif' WRITE(inum,*) ' total changes among all mpi tasks:' WRITE(inum,*) ' mpiwe mpiea mpiso mpino mpisw mpise mpinw mpine' WRITE(inum,'(a,8i6)') ' Send: ', ichanged(jh*16-15:jh*16-8) WRITE(inum,'(a,8i6)') ' Recv: ', ichanged(jh*16 -7:jh*16 ) END DO ENDIF ! CALL init_ioipsl ! Prepare NetCDF output file (if necessary) ! IF (llwrtlay) CLOSE(inum) ! DEALLOCATE(iin, ijn, iimppt, ijmppt, ijpi, ijpj, ipproc, inei, llnei, impi, llisOce) ! END SUBROUTINE mpp_init #endif SUBROUTINE mpp_basesplit( kiglo, kjglo, khls, knbi, knbj, kimax, kjmax, kimppt, kjmppt, klci, klcj) !!---------------------------------------------------------------------- !! *** ROUTINE mpp_basesplit *** !! !! ** Purpose : Lay out the global domain over processors. !! !! ** Method : Global domain is distributed in smaller local domains. !! !! ** Action : - set for all knbi*knbj domains: !! kimppt : longitudinal index !! kjmppt : latitudinal index !! klci : first dimension !! klcj : second dimension !!---------------------------------------------------------------------- INTEGER, INTENT(in ) :: kiglo, kjglo INTEGER, INTENT(in ) :: khls INTEGER, INTENT(in ) :: knbi, knbj INTEGER, INTENT( out) :: kimax, kjmax INTEGER, DIMENSION(knbi,knbj), OPTIONAL, INTENT( out) :: kimppt, kjmppt INTEGER, DIMENSION(knbi,knbj), OPTIONAL, INTENT( out) :: klci, klcj ! INTEGER :: ji, jj INTEGER :: i2hls INTEGER :: iresti, irestj, irm, ijpjmin !!---------------------------------------------------------------------- i2hls = 2*khls ! #if defined key_nemocice_decomp kimax = ( nx_global+2-i2hls + (knbi-1) ) / knbi + i2hls ! first dim. kjmax = ( ny_global+2-i2hls + (knbj-1) ) / knbj + i2hls ! second dim. #else kimax = ( kiglo - i2hls + (knbi-1) ) / knbi + i2hls ! first dim. kjmax = ( kjglo - i2hls + (knbj-1) ) / knbj + i2hls ! second dim. #endif IF( .NOT. PRESENT(kimppt) ) RETURN ! ! 1. Dimension arrays for subdomains ! ----------------------------------- ! Computation of local domain sizes klci() klcj() ! These dimensions depend on global sizes knbi,knbj and kiglo,kjglo ! The subdomains are squares lesser than or equal to the global ! dimensions divided by the number of processors minus the overlap array. ! iresti = 1 + MOD( kiglo - i2hls - 1 , knbi ) irestj = 1 + MOD( kjglo - i2hls - 1 , knbj ) ! ! Need to use kimax and kjmax here since jpi and jpj not yet defined #if defined key_nemocice_decomp ! Change padding to be consistent with CICE klci(1:knbi-1,: ) = kimax klci( knbi ,: ) = kiglo - (knbi - 1) * (kimax - i2hls) klcj(: ,1:knbj-1) = kjmax klcj(: , knbj ) = kjglo - (knbj - 1) * (kjmax - i2hls) #else klci(1:iresti ,:) = kimax klci(iresti+1:knbi ,:) = kimax-1 IF( MINVAL(klci) < 2*i2hls ) THEN WRITE(ctmp1,*) ' mpp_basesplit: minimum value of jpi must be >= ', 2*i2hls WRITE(ctmp2,*) ' We have ', MINVAL(klci) CALL ctl_stop( 'STOP', ctmp1, ctmp2 ) ENDIF IF( l_NFold ) THEN ! minimize the size of the last row to compensate for the north pole folding coast IF( c_NFtype == 'T' ) ijpjmin = 2+3*khls ! V and F folding must be outside of southern halos IF( c_NFtype == 'F' ) ijpjmin = 1+3*khls ! V and F folding must be outside of southern halos irm = knbj - irestj ! total number of lines to be removed klcj(:,knbj) = MAX( ijpjmin, kjmax-irm ) ! we must have jpj >= ijpjmin in the last row irm = irm - ( kjmax - klcj(1,knbj) ) ! remaining number of lines to remove irestj = knbj - 1 - irm klcj(:, irestj+1:knbj-1) = kjmax-1 ELSE klcj(:, irestj+1:knbj ) = kjmax-1 ENDIF klcj(:,1:irestj) = kjmax IF( MINVAL(klcj) < 2*i2hls ) THEN WRITE(ctmp1,*) ' mpp_basesplit: minimum value of jpj must be >= ', 2*i2hls WRITE(ctmp2,*) ' We have ', MINVAL(klcj) CALL ctl_stop( 'STOP', ctmp1, ctmp2 ) ENDIF #endif ! 2. Index arrays for subdomains ! ------------------------------- kimppt(:,:) = 1 kjmppt(:,:) = 1 ! IF( knbi > 1 ) THEN DO jj = 1, knbj DO ji = 2, knbi kimppt(ji,jj) = kimppt(ji-1,jj) + klci(ji-1,jj) - i2hls END DO END DO ENDIF ! IF( knbj > 1 )THEN DO jj = 2, knbj DO ji = 1, knbi kjmppt(ji,jj) = kjmppt(ji,jj-1) + klcj(ji,jj-1) - i2hls END DO END DO ENDIF END SUBROUTINE mpp_basesplit SUBROUTINE bestpartition( knbij, knbi, knbj, knbcnt, ldlist ) !!---------------------------------------------------------------------- !! *** ROUTINE bestpartition *** !! !! ** Purpose : !! !! ** Method : !!---------------------------------------------------------------------- INTEGER, INTENT(in ) :: knbij ! total number of subdomains (knbi*knbj) INTEGER, OPTIONAL, INTENT( out) :: knbi, knbj ! number if subdomains along i and j (knbi and knbj) INTEGER, OPTIONAL, INTENT( out) :: knbcnt ! number of land subdomains LOGICAL, OPTIONAL, INTENT(in ) :: ldlist ! .true.: print the list the best domain decompositions (with land) ! INTEGER :: ji, jj, ii, iitarget INTEGER :: iszitst, iszjtst INTEGER :: isziref, iszjref INTEGER :: iszimin, iszjmin INTEGER :: inbij, iszij INTEGER :: inbimax, inbjmax, inbijmax, inbijold INTEGER :: isz0, isz1 INTEGER, DIMENSION( :), ALLOCATABLE :: indexok INTEGER, DIMENSION( :), ALLOCATABLE :: inbi0, inbj0, inbij0 ! number of subdomains along i,j INTEGER, DIMENSION( :), ALLOCATABLE :: iszi0, iszj0, iszij0 ! max size of the subdomains along i,j INTEGER, DIMENSION( :), ALLOCATABLE :: inbi1, inbj1, inbij1 ! number of subdomains along i,j INTEGER, DIMENSION( :), ALLOCATABLE :: iszi1, iszj1, iszij1 ! max size of the subdomains along i,j LOGICAL :: llist LOGICAL, DIMENSION(:,:), ALLOCATABLE :: llmsk2d ! max size of the subdomains along i,j LOGICAL, DIMENSION(:,:), ALLOCATABLE :: llisOce ! - - REAL(wp):: zpropland !!---------------------------------------------------------------------- ! llist = .FALSE. IF( PRESENT(ldlist) ) llist = ldlist CALL mpp_init_landprop( zpropland ) ! get the proportion of land point over the gloal domain inbij = NINT( REAL(knbij, wp) / ( 1.0 - zpropland ) ) ! define the largest possible value for jpni*jpnj ! IF( llist ) THEN ; inbijmax = inbij*2 ELSE ; inbijmax = inbij ENDIF ! ALLOCATE(inbi0(inbijmax),inbj0(inbijmax),iszi0(inbijmax),iszj0(inbijmax)) ! inbimax = 0 inbjmax = 0 isziref = jpiglo*jpjglo+1 ! define a value that is larger than the largest possible iszjref = jpiglo*jpjglo+1 ! iszimin = 4*nn_hls ! minimum size of the MPI subdomain so halos are always adressing neighbor inner domain iszjmin = 4*nn_hls IF( c_NFtype == 'T' ) iszjmin = MAX(iszjmin, 2+3*nn_hls) ! V and F folding must be outside of southern halos IF( c_NFtype == 'F' ) iszjmin = MAX(iszjmin, 1+3*nn_hls) ! V and F folding must be outside of southern halos ! ! get the list of knbi that gives a smaller jpimax than knbi-1 ! get the list of knbj that gives a smaller jpjmax than knbj-1 DO ji = 1, inbijmax #if defined key_nemocice_decomp iszitst = ( nx_global+2-2*nn_hls + (ji-1) ) / ji + 2*nn_hls ! first dim. #else iszitst = ( Ni0glo + (ji-1) ) / ji + 2*nn_hls ! max subdomain i-size #endif IF( iszitst < isziref .AND. iszitst >= iszimin ) THEN isziref = iszitst inbimax = inbimax + 1 inbi0(inbimax) = ji iszi0(inbimax) = isziref ENDIF #if defined key_nemocice_decomp iszjtst = ( ny_global+2-2*nn_hls + (ji-1) ) / ji + 2*nn_hls ! first dim. #else iszjtst = ( Nj0glo + (ji-1) ) / ji + 2*nn_hls ! max subdomain j-size #endif IF( iszjtst < iszjref .AND. iszjtst >= iszjmin ) THEN iszjref = iszjtst inbjmax = inbjmax + 1 inbj0(inbjmax) = ji iszj0(inbjmax) = iszjref ENDIF END DO ! combine these 2 lists to get all possible knbi*knbj < inbijmax ALLOCATE( llmsk2d(inbimax,inbjmax) ) DO jj = 1, inbjmax DO ji = 1, inbimax IF ( inbi0(ji) * inbj0(jj) <= inbijmax ) THEN ; llmsk2d(ji,jj) = .TRUE. ELSE ; llmsk2d(ji,jj) = .FALSE. ENDIF END DO END DO isz1 = COUNT(llmsk2d) ALLOCATE( inbi1(isz1), inbj1(isz1), iszi1(isz1), iszj1(isz1) ) ii = 0 DO jj = 1, inbjmax DO ji = 1, inbimax IF( llmsk2d(ji,jj) .EQV. .TRUE. ) THEN ii = ii + 1 inbi1(ii) = inbi0(ji) inbj1(ii) = inbj0(jj) iszi1(ii) = iszi0(ji) iszj1(ii) = iszj0(jj) ENDIF END DO END DO DEALLOCATE( inbi0, inbj0, iszi0, iszj0 ) DEALLOCATE( llmsk2d ) ALLOCATE( inbij1(isz1), iszij1(isz1) ) inbij1(:) = inbi1(:) * inbj1(:) iszij1(:) = iszi1(:) * iszj1(:) ! if there is no land and no print IF( .NOT. llist .AND. numbot == -1 .AND. numbdy == -1 ) THEN ! get the smaller partition which gives the smallest subdomain size ii = MINLOC(inbij1, mask = iszij1 == MINVAL(iszij1), dim = 1) knbi = inbi1(ii) knbj = inbj1(ii) IF(PRESENT(knbcnt)) knbcnt = 0 DEALLOCATE( inbi1, inbj1, inbij1, iszi1, iszj1, iszij1 ) RETURN ENDIF ! extract only the partitions which reduce the subdomain size in comparison with smaller partitions ALLOCATE( indexok(isz1) ) ! to store indices of the best partitions isz0 = 0 ! number of best partitions inbij = 1 ! start with the min value of inbij1 => 1 iszij = jpiglo*jpjglo+1 ! default: larger than global domain DO WHILE( inbij <= inbijmax ) ! if we did not reach the max of inbij1 ii = MINLOC(iszij1, mask = inbij1 == inbij, dim = 1) ! warning: send back the first occurence if multiple results IF ( iszij1(ii) < iszij ) THEN ii = MINLOC( iszi1+iszj1, mask = iszij1 == iszij1(ii) .AND. inbij1 == inbij, dim = 1) ! select the smaller perimeter if multiple min isz0 = isz0 + 1 indexok(isz0) = ii iszij = iszij1(ii) ENDIF inbij = MINVAL(inbij1, mask = inbij1 > inbij) ! warning: return largest integer value if mask = .false. everywhere END DO DEALLOCATE( inbij1, iszij1 ) ! keep only the best partitions (sorted by increasing order of subdomains number and decreassing subdomain size) ALLOCATE( inbi0(isz0), inbj0(isz0), iszi0(isz0), iszj0(isz0) ) DO ji = 1, isz0 ii = indexok(ji) inbi0(ji) = inbi1(ii) inbj0(ji) = inbj1(ii) iszi0(ji) = iszi1(ii) iszj0(ji) = iszj1(ii) END DO DEALLOCATE( indexok, inbi1, inbj1, iszi1, iszj1 ) IF( llist ) THEN IF(lwp) THEN WRITE(numout,*) WRITE(numout,*) ' For your information:' WRITE(numout,*) ' list of the best partitions including land supression' WRITE(numout,*) ' -----------------------------------------------------' WRITE(numout,*) ENDIF ji = isz0 ! initialization with the largest value ALLOCATE( llisOce(inbi0(ji), inbj0(ji)) ) CALL mpp_is_ocean( llisOce ) ! Warning: must be call by all cores (call mpp_sum) inbijold = COUNT(llisOce) DEALLOCATE( llisOce ) DO ji =isz0-1,1,-1 ALLOCATE( llisOce(inbi0(ji), inbj0(ji)) ) CALL mpp_is_ocean( llisOce ) ! Warning: must be call by all cores (call mpp_sum) inbij = COUNT(llisOce) DEALLOCATE( llisOce ) IF(lwp .AND. inbij < inbijold) THEN WRITE(numout,'(a, i6, a, i6, a, f4.1, a, i9, a, i6, a, i6, a)') & & 'nb_cores oce: ', inbij, ', land domains excluded: ', inbi0(ji)*inbj0(ji) - inbij, & & ' (', REAL(inbi0(ji)*inbj0(ji) - inbij,wp) / REAL(inbi0(ji)*inbj0(ji),wp) *100., & & '%), largest oce domain: ', iszi0(ji)*iszj0(ji), ' ( ', iszi0(ji),' x ', iszj0(ji), ' )' inbijold = inbij ENDIF END DO DEALLOCATE( inbi0, inbj0, iszi0, iszj0 ) IF(lwp) THEN WRITE(numout,*) WRITE(numout,*) ' -----------------------------------------------------------' ENDIF CALL mppsync CALL mppstop( ld_abort = .TRUE. ) ENDIF DEALLOCATE( iszi0, iszj0 ) inbij = inbijmax + 1 ! default: larger than possible ii = isz0+1 ! start from the end of the list (smaller subdomains) DO WHILE( inbij > knbij ) ! while the number of ocean subdomains exceed the number of procs ii = ii -1 ALLOCATE( llisOce(inbi0(ii), inbj0(ii)) ) CALL mpp_is_ocean( llisOce ) ! must be done by all core inbij = COUNT(llisOce) DEALLOCATE( llisOce ) END DO knbi = inbi0(ii) knbj = inbj0(ii) IF(PRESENT(knbcnt)) knbcnt = knbi * knbj - inbij DEALLOCATE( inbi0, inbj0 ) ! END SUBROUTINE bestpartition SUBROUTINE mpp_init_landprop( propland ) !!---------------------------------------------------------------------- !! *** ROUTINE mpp_init_landprop *** !! !! ** Purpose : the the proportion of land points in the surface land-sea mask !! !! ** Method : read iproc strips (of length Ni0glo) of the land-sea mask !!---------------------------------------------------------------------- REAL(wp), INTENT( out) :: propland ! proportion of land points in the global domain (between 0 and 1) ! INTEGER, DIMENSION(jpni*jpnj) :: kusedom_1d INTEGER :: inboce, iarea INTEGER :: iproc, idiv, ijsz INTEGER :: ijstr LOGICAL, ALLOCATABLE, DIMENSION(:,:) :: lloce !!---------------------------------------------------------------------- ! do nothing if there is no land-sea mask IF( numbot == -1 .and. numbdy == -1 ) THEN propland = 0. RETURN ENDIF ! number of processes reading the bathymetry file iproc = MINVAL( (/mppsize, Nj0glo/2, 100/) ) ! read a least 2 lines, no more that 100 processes reading at the same time ! we want to read iproc strips of the land-sea mask. -> pick up iproc processes every idiv processes starting at 1 IF( iproc == 1 ) THEN ; idiv = mppsize ELSE ; idiv = ( mppsize - 1 ) / ( iproc - 1 ) ENDIF iarea = (narea-1)/idiv ! involed process number (starting counting at 0) IF( MOD( narea-1, idiv ) == 0 .AND. iarea < iproc ) THEN ! beware idiv can be = to 1 ! ijsz = Nj0glo / iproc ! width of the stripe to read IF( iarea < MOD(Nj0glo,iproc) ) ijsz = ijsz + 1 ijstr = iarea*(Nj0glo/iproc) + MIN(iarea, MOD(Nj0glo,iproc)) + 1 ! starting j position of the reading ! ALLOCATE( lloce(Ni0glo, ijsz) ) ! allocate the strip CALL read_mask( 1, ijstr, Ni0glo, ijsz, lloce ) inboce = COUNT(lloce) ! number of ocean point in the stripe DEALLOCATE(lloce) ! ELSE inboce = 0 ENDIF CALL mpp_sum( 'mppini', inboce ) ! total number of ocean points over the global domain ! propland = REAL( Ni0glo*Nj0glo - inboce, wp ) / REAL( Ni0glo*Nj0glo, wp ) ! END SUBROUTINE mpp_init_landprop SUBROUTINE mpp_is_ocean( ldIsOce ) !!---------------------------------------------------------------------- !! *** ROUTINE mpp_is_ocean *** !! !! ** Purpose : Check for a mpi domain decomposition inbi x inbj which !! subdomains, including 1 halo (even if nn_hls>1), contain !! at least 1 ocean point. !! We must indeed ensure that each subdomain that is a neighbour !! of a land subdomain, has only land points on its boundary !! (inside the inner subdomain) with the land subdomain. !! This is needed to get the proper bondary conditions on !! a subdomain with a closed boundary. !! !! ** Method : read inbj strips (of length Ni0glo) of the land-sea mask !!---------------------------------------------------------------------- LOGICAL, DIMENSION(:,:), INTENT( out) :: ldIsOce ! .true. if a sub domain constains 1 ocean point ! INTEGER :: idiv, iimax, ijmax, iarea INTEGER :: inbi, inbj, inx, iny, inry, isty INTEGER :: ji, jn INTEGER, ALLOCATABLE, DIMENSION(:,:) :: inboce ! number oce oce pint in each mpi subdomain INTEGER, ALLOCATABLE, DIMENSION(: ) :: inboce_1d INTEGER, ALLOCATABLE, DIMENSION(:,:) :: iimppt, ijpi INTEGER, ALLOCATABLE, DIMENSION(:,:) :: ijmppt, ijpj LOGICAL, ALLOCATABLE, DIMENSION(:,:) :: lloce ! lloce(i,j) = .true. if the point (i,j) is ocean !!---------------------------------------------------------------------- ! do nothing if there is no land-sea mask IF( numbot == -1 .AND. numbdy == -1 ) THEN ldIsOce(:,:) = .TRUE. RETURN ENDIF ! inbi = SIZE( ldIsOce, dim = 1 ) inbj = SIZE( ldIsOce, dim = 2 ) ! ! we want to read inbj strips of the land-sea mask. -> pick up inbj processes every idiv processes starting at 1 IF ( inbj == 1 ) THEN ; idiv = mppsize ELSE IF ( mppsize < inbj ) THEN ; idiv = 1 ELSE ; idiv = ( mppsize - 1 ) / ( inbj - 1 ) ENDIF ! ALLOCATE( inboce(inbi,inbj), inboce_1d(inbi*inbj) ) inboce(:,:) = 0 ! default no ocean point found ! DO jn = 0, (inbj-1)/mppsize ! if mppsize < inbj : more strips than mpi processes (because of potential land domains) ! iarea = (narea-1)/idiv + jn * mppsize + 1 ! involed process number (starting counting at 1) IF( MOD( narea-1, idiv ) == 0 .AND. iarea <= inbj ) THEN ! beware idiv can be = to 1 ! ALLOCATE( iimppt(inbi,inbj), ijmppt(inbi,inbj), ijpi(inbi,inbj), ijpj(inbi,inbj) ) CALL mpp_basesplit( Ni0glo, Nj0glo, 0, inbi, inbj, iimax, ijmax, iimppt, ijmppt, ijpi, ijpj ) ! inx = Ni0glo + 2 ; iny = ijpj(1,iarea) + 2 ! strip size + 1 halo on each direction (even if nn_hls>1) ALLOCATE( lloce(inx, iny) ) ! allocate the strip inry = iny - COUNT( (/ iarea == 1, iarea == inbj /) ) ! number of point to read in y-direction isty = 1 + COUNT( (/ iarea == 1 /) ) ! read from the first or the second line? CALL read_mask( 1, ijmppt(1,iarea) - 2 + isty, Ni0glo, inry, lloce(2:inx-1, isty:inry+isty-1) ) ! read the strip ! IF( iarea == 1 ) THEN ! the first line was not read IF( l_Jperio ) THEN ! north-south periodocity CALL read_mask( 1, Nj0glo, Ni0glo, 1, lloce(2:inx-1, 1) ) ! read the last line -> first line of lloce ELSE lloce(2:inx-1, 1) = .FALSE. ! closed boundary ENDIF ENDIF IF( iarea == inbj ) THEN ! the last line was not read IF( l_Jperio ) THEN ! north-south periodocity CALL read_mask( 1, 1, Ni0glo, 1, lloce(2:inx-1,iny) ) ! read the first line -> last line of lloce ELSEIF( c_NFtype == 'T' ) THEN ! north-pole folding T-pivot, T-point lloce(2,iny) = lloce(2,iny-2) ! here we have 1 halo (even if nn_hls>1) DO ji = 3,inx-1 lloce(ji,iny ) = lloce(inx-ji+2,iny-2) ! ok, we have at least 3 lines END DO DO ji = inx/2+2,inx-1 lloce(ji,iny-1) = lloce(inx-ji+2,iny-1) END DO ELSEIF( c_NFtype == 'F' ) THEN ! north-pole folding F-pivot, T-point, 1 halo lloce(inx/2+1,iny-1) = lloce(inx/2,iny-1) ! here we have 1 halo (even if nn_hls>1) lloce(inx -1,iny-1) = lloce(2 ,iny-1) DO ji = 2,inx-1 lloce(ji,iny) = lloce(inx-ji+1,iny-1) END DO ELSE ! closed boundary lloce(2:inx-1,iny) = .FALSE. ENDIF ENDIF ! ! first and last column were not read IF( l_Iperio ) THEN lloce(1,:) = lloce(inx-1,:) ; lloce(inx,:) = lloce(2,:) ! east-west periodocity ELSE lloce(1,:) = .FALSE. ; lloce(inx,:) = .FALSE. ! closed boundary ENDIF ! DO ji = 1, inbi inboce(ji,iarea) = COUNT( lloce(iimppt(ji,1):iimppt(ji,1)+ijpi(ji,1)+1,:) ) ! lloce as 2 points more than Ni0glo END DO ! DEALLOCATE(lloce) DEALLOCATE(iimppt, ijmppt, ijpi, ijpj) ! ENDIF END DO inboce_1d = RESHAPE(inboce, (/ inbi*inbj /)) CALL mpp_sum( 'mppini', inboce_1d ) inboce = RESHAPE(inboce_1d, (/inbi, inbj/)) ldIsOce(:,:) = inboce(:,:) /= 0 DEALLOCATE(inboce, inboce_1d) ! END SUBROUTINE mpp_is_ocean SUBROUTINE read_mask( kistr, kjstr, kicnt, kjcnt, ldoce ) !!---------------------------------------------------------------------- !! *** ROUTINE read_mask *** !! !! ** Purpose : Read relevant bathymetric information in order to !! provide a land/sea mask used for the elimination !! of land domains, in an mpp computation. !! !! ** Method : read stipe of size (Ni0glo,...) !!---------------------------------------------------------------------- INTEGER , INTENT(in ) :: kistr, kjstr ! starting i and j position of the reading INTEGER , INTENT(in ) :: kicnt, kjcnt ! number of points to read in i and j directions LOGICAL, DIMENSION(kicnt,kjcnt), INTENT( out) :: ldoce ! ldoce(i,j) = .true. if the point (i,j) is ocean ! INTEGER :: inumsave ! local logical unit REAL(wp), DIMENSION(kicnt,kjcnt) :: zbot, zbdy !!---------------------------------------------------------------------- ! inumsave = numout ; numout = numnul ! redirect all print to /dev/null ! IF( numbot /= -1 ) THEN CALL iom_get( numbot, jpdom_unknown, 'bottom_level', zbot, kstart = (/kistr,kjstr/), kcount = (/kicnt, kjcnt/) ) ELSE zbot(:,:) = 1._wp ! put a non-null value ENDIF ! IF( numbdy /= -1 ) THEN ! Adjust with bdy_msk if it exists CALL iom_get ( numbdy, jpdom_unknown, 'bdy_msk', zbdy, kstart = (/kistr,kjstr/), kcount = (/kicnt, kjcnt/) ) zbot(:,:) = zbot(:,:) * zbdy(:,:) ENDIF ! ldoce(:,:) = NINT(zbot(:,:)) > 0 numout = inumsave ! END SUBROUTINE read_mask SUBROUTINE mpp_getnum( ldIsOce, kproc, kipos, kjpos ) !!---------------------------------------------------------------------- !! *** ROUTINE mpp_getnum *** !! !! ** Purpose : give a number to each MPI subdomains (starting at 0) !! !! ** Method : start from bottom left. First skip land subdomain, and finally use them if needed !!---------------------------------------------------------------------- LOGICAL, DIMENSION(:,:), INTENT(in ) :: ldIsOce ! F if land process INTEGER, DIMENSION(:,:), INTENT( out) :: kproc ! subdomain number (-1 if not existing, starting at 0) INTEGER, DIMENSION( :), INTENT( out) :: kipos ! i-position of the subdomain (from 1 to jpni) INTEGER, DIMENSION( :), INTENT( out) :: kjpos ! j-position of the subdomain (from 1 to jpnj) ! INTEGER :: ii, ij, jarea, iarea0 INTEGER :: icont, i2add , ini, inj, inij !!---------------------------------------------------------------------- ! ini = SIZE(ldIsOce, dim = 1) inj = SIZE(ldIsOce, dim = 2) inij = SIZE(kipos) ! ! specify which subdomains are oce subdomains; other are land subdomains kproc(:,:) = -1 icont = -1 DO jarea = 1, ini*inj iarea0 = jarea - 1 ii = 1 + MOD(iarea0,ini) ij = 1 + iarea0/ini IF( ldIsOce(ii,ij) ) THEN icont = icont + 1 kproc(ii,ij) = icont kipos(icont+1) = ii kjpos(icont+1) = ij ENDIF END DO ! if needed add some land subdomains to reach inij active subdomains i2add = inij - COUNT( ldIsOce ) DO jarea = 1, ini*inj iarea0 = jarea - 1 ii = 1 + MOD(iarea0,ini) ij = 1 + iarea0/ini IF( .NOT. ldIsOce(ii,ij) .AND. i2add > 0 ) THEN icont = icont + 1 kproc(ii,ij) = icont kipos(icont+1) = ii kjpos(icont+1) = ij i2add = i2add - 1 ENDIF END DO ! END SUBROUTINE mpp_getnum SUBROUTINE init_excl_landpt !!---------------------------------------------------------------------- !! *** ROUTINE *** !! !! ** Purpose : exclude exchanges which contain only land points !! !! ** Method : if a send or receive buffer constains only land point we !! flag off the corresponding communication !! Warning: this selection depend on the halo size -> loop on halo size !! !!---------------------------------------------------------------------- INTEGER :: inumsave INTEGER :: jh INTEGER :: ipi, ipj INTEGER :: iiwe, iiea, iist, iisz INTEGER :: ijso, ijno, ijst, ijsz LOGICAL :: llsave REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zmsk LOGICAL , DIMENSION(Ni_0,Nj_0,1) :: lloce !!---------------------------------------------------------------------- ! ! read the land-sea mask on the inner domain CALL read_mask( nimpp, njmpp, Ni_0, Nj_0, lloce(:,:,1) ) ! ! Here we look only at communications excluding the NP folding. ! As lbcnfd not validated if halo size /= nn_hls, we switch if off temporary... llsave = l_IdoNFold l_IdoNFold = .FALSE. ! DO jh = 1, n_hlsmax ! different halo size ! ipi = Ni_0 + 2*jh ! local domain size ipj = Nj_0 + 2*jh ! ALLOCATE( zmsk(ipi,ipj) ) 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 CALL lbc_lnk('mppini', zmsk, 'T', 1._wp, khls = jh) ! fill halos ! iiwe = jh ; iiea = Ni_0 ! bottom-left corfer - 1 of the sent data ijso = jh ; ijno = Nj_0 IF( nn_comm == 1 ) THEN iist = 0 ; iisz = ipi ijst = 0 ; ijsz = ipj ELSE iist = jh ; iisz = Ni_0 ijst = jh ; ijsz = Nj_0 ENDIF IF( nn_comm == 1 ) THEN ! SM: NOT WORKING FOR NEIGHBOURHOOD COLLECTIVE COMMUNICATIONS, I DON'T KNOW WHY... ! do not send if we send only land points IF( NINT(SUM( zmsk(iiwe+1:iiwe+jh ,ijst+1:ijst+ijsz) )) == 0 ) mpiSnei(jh,jpwe) = -1 IF( NINT(SUM( zmsk(iiea+1:iiea+jh ,ijst+1:ijst+ijsz) )) == 0 ) mpiSnei(jh,jpea) = -1 IF( NINT(SUM( zmsk(iist+1:iist+iisz,ijso+1:ijso+jh ) )) == 0 ) mpiSnei(jh,jpso) = -1 IF( NINT(SUM( zmsk(iist+1:iist+iisz,ijno+1:ijno+jh ) )) == 0 ) mpiSnei(jh,jpno) = -1 IF( NINT(SUM( zmsk(iiwe+1:iiwe+jh ,ijso+1:ijso+jh ) )) == 0 ) mpiSnei(jh,jpsw) = -1 IF( NINT(SUM( zmsk(iiea+1:iiea+jh ,ijso+1:ijso+jh ) )) == 0 ) mpiSnei(jh,jpse) = -1 IF( NINT(SUM( zmsk(iiwe+1:iiwe+jh ,ijno+1:ijno+jh ) )) == 0 ) mpiSnei(jh,jpnw) = -1 IF( NINT(SUM( zmsk(iiea+1:iiea+jh ,ijno+1:ijno+jh ) )) == 0 ) mpiSnei(jh,jpne) = -1 ! iiwe = iiwe-jh ; iiea = iiea+jh ! bottom-left corfer - 1 of the received data ijso = ijso-jh ; ijno = ijno+jh ! do not send if we send only land points IF( NINT(SUM( zmsk(iiwe+1:iiwe+jh ,ijst+1:ijst+ijsz) )) == 0 ) mpiRnei(jh,jpwe) = -1 IF( NINT(SUM( zmsk(iiea+1:iiea+jh ,ijst+1:ijst+ijsz) )) == 0 ) mpiRnei(jh,jpea) = -1 IF( NINT(SUM( zmsk(iist+1:iist+iisz,ijso+1:ijso+jh ) )) == 0 ) mpiRnei(jh,jpso) = -1 IF( NINT(SUM( zmsk(iist+1:iist+iisz,ijno+1:ijno+jh ) )) == 0 ) mpiRnei(jh,jpno) = -1 IF( NINT(SUM( zmsk(iiwe+1:iiwe+jh ,ijso+1:ijso+jh ) )) == 0 ) mpiRnei(jh,jpsw) = -1 IF( NINT(SUM( zmsk(iiea+1:iiea+jh ,ijso+1:ijso+jh ) )) == 0 ) mpiRnei(jh,jpse) = -1 IF( NINT(SUM( zmsk(iiwe+1:iiwe+jh ,ijno+1:ijno+jh ) )) == 0 ) mpiRnei(jh,jpnw) = -1 IF( NINT(SUM( zmsk(iiea+1:iiea+jh ,ijno+1:ijno+jh ) )) == 0 ) mpiRnei(jh,jpne) = -1 ENDIF ! ! Specific (and rare) problem in corner treatment because we do 1st West-East comm, next South-North comm IF( nn_comm == 1 ) THEN IF( mpiSnei(jh,jpwe) > -1 ) mpiSnei(jh, (/jpsw,jpnw/) ) = -1 ! SW and NW corners already sent through West nei IF( mpiSnei(jh,jpea) > -1 ) mpiSnei(jh, (/jpse,jpne/) ) = -1 ! SE and NE corners already sent through East nei IF( mpiRnei(jh,jpso) > -1 ) mpiRnei(jh, (/jpsw,jpse/) ) = -1 ! SW and SE corners will be received through South nei IF( mpiRnei(jh,jpno) > -1 ) mpiRnei(jh, (/jpnw,jpne/) ) = -1 ! NW and NE corners will be received through North nei ENDIF ! DEALLOCATE( zmsk ) ! CALL mpp_ini_nc(jh) ! Initialize/Update communicator for neighbourhood collective communications ! END DO l_IdoNFold = llsave END SUBROUTINE init_excl_landpt SUBROUTINE init_ioipsl !!---------------------------------------------------------------------- !! *** ROUTINE init_ioipsl *** !! !! ** Purpose : !! !! ** Method : !! !! History : !! 9.0 ! 04-03 (G. Madec ) MPP-IOIPSL !! " " ! 08-12 (A. Coward) addition in case of jpni*jpnj < jpnij !!---------------------------------------------------------------------- INTEGER, DIMENSION(2) :: iglo, iloc, iabsf, iabsl, ihals, ihale, idid !!---------------------------------------------------------------------- ! The domain is split only horizontally along i- or/and j- direction ! So we need at the most only 1D arrays with 2 elements. ! Set idompar values equivalent to the jpdom_local_noextra definition ! used in IOM. This works even if jpnij .ne. jpni*jpnj. iglo( :) = (/ Ni0glo, Nj0glo /) iloc( :) = (/ Ni_0 , Nj_0 /) iabsf(:) = (/ Nis0 , Njs0 /) + (/ nimpp, njmpp /) - 1 - nn_hls ! corresponds to mig0(Nis0) but mig0 is not yet defined! iabsl(:) = iabsf(:) + iloc(:) - 1 ihals(:) = (/ 0 , 0 /) ihale(:) = (/ 0 , 0 /) idid( :) = (/ 1 , 2 /) IF(lwp) THEN WRITE(numout,*) WRITE(numout,*) 'mpp init_ioipsl : iloc = ', iloc WRITE(numout,*) '~~~~~~~~~~~~~~~ iabsf = ', iabsf WRITE(numout,*) ' ihals = ', ihals WRITE(numout,*) ' ihale = ', ihale ENDIF ! CALL flio_dom_set ( jpnij, narea-1, idid, iglo, iloc, iabsf, iabsl, ihals, ihale, 'BOX', nidom) ! END SUBROUTINE init_ioipsl SUBROUTINE init_nfdcom !!---------------------------------------------------------------------- !! *** ROUTINE init_nfdcom *** !! ** Purpose : Setup for north fold exchanges with explicit !! point-to-point messaging !! !! ** Method : Initialization of the northern neighbours lists. !!---------------------------------------------------------------------- !! 1.0 ! 2011-10 (A. C. Coward, NOCS & J. Donners, PRACE) !! 2.0 ! 2013-06 Setup avoiding MPI communication (I. Epicoco, S. Mocavero, CMCC) !!---------------------------------------------------------------------- INTEGER :: sxM, dxM, sxT, dxT, jn !!---------------------------------------------------------------------- ! !sxM is the first point (in the global domain) needed to compute the north-fold for the current process sxM = jpiglo - nimpp - jpi + 1 !dxM is the last point (in the global domain) needed to compute the north-fold for the current process dxM = jpiglo - nimpp + 2 ! ! loop over the other north-fold processes to find the processes ! managing the points belonging to the sxT-dxT range ! nsndto = 0 DO jn = 1, jpni ! sxT = nfimpp(jn) ! sxT = 1st point (in the global domain) of the jn process dxT = nfimpp(jn) + nfjpi(jn) - 1 ! dxT = last point (in the global domain) of the jn process ! IF ( sxT < sxM .AND. sxM < dxT ) THEN nsndto = nsndto + 1 isendto(nsndto) = jn ELSEIF( sxM <= sxT .AND. dxM >= dxT ) THEN nsndto = nsndto + 1 isendto(nsndto) = jn ELSEIF( dxM < dxT .AND. sxT < dxM ) THEN nsndto = nsndto + 1 isendto(nsndto) = jn ENDIF ! END DO ! END SUBROUTINE init_nfdcom SUBROUTINE init_doloop !!---------------------------------------------------------------------- !! *** ROUTINE init_doloop *** !! !! ** Purpose : set the starting/ending indices of DO-loop !! These indices are used in do_loop_substitute.h90 !!---------------------------------------------------------------------- ! Nis0 = 1+nn_hls Njs0 = 1+nn_hls Nie0 = jpi-nn_hls Nje0 = jpj-nn_hls ! Ni_0 = Nie0 - Nis0 + 1 Nj_0 = Nje0 - Njs0 + 1 ! ! old indices to be removed... jpim1 = jpi-1 ! inner domain indices jpjm1 = jpj-1 ! " " jpkm1 = jpk-1 ! " " ! END SUBROUTINE init_doloop !!====================================================================== END MODULE mppini