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.
lbc_lnk_neicoll_generic.h90 in NEMO/trunk/src/OCE/LBC – NEMO

source: NEMO/trunk/src/OCE/LBC/lbc_lnk_neicoll_generic.h90 @ 15725

Last change on this file since 15725 was 15354, checked in by smasson, 3 years ago

trunk: BDY compliant with corner communications, see #2731

File size: 15.0 KB
RevLine 
[14021]1
[14366]2   SUBROUTINE lbc_lnk_neicoll_/**/PRECISION( cdname, ptab, cd_nat, psgn, kfld, kfillmode, pfillval, khls, lsend, lrecv, ld4only )
[14021]3      CHARACTER(len=*)              , INTENT(in   ) ::   cdname      ! name of the calling subroutine
[14349]4      TYPE(PTR_4d_/**/PRECISION),  DIMENSION(:), INTENT(inout) ::   ptab        ! pointer of arrays on which apply the b.c.
[14338]5      CHARACTER(len=1), DIMENSION(:), INTENT(in   ) ::   cd_nat      ! nature of array grid-points
6      REAL(PRECISION),  DIMENSION(:), INTENT(in   ) ::   psgn        ! sign used across the north fold boundary
7      INTEGER                       , INTENT(in   ) ::   kfld        ! number of pt3d arrays
[14021]8      INTEGER ,             OPTIONAL, INTENT(in   ) ::   kfillmode   ! filling method for halo over land (default = constant)
[14314]9      REAL(PRECISION),      OPTIONAL, INTENT(in   ) ::   pfillval    ! background value (used at closed boundaries)
[14363]10      INTEGER ,             OPTIONAL, INTENT(in   ) ::   khls        ! halo size, default = nn_hls
[15354]11      LOGICAL, DIMENSION(8),OPTIONAL, INTENT(in   ) ::   lsend, lrecv  ! communication with other 4 proc
[14366]12      LOGICAL,              OPTIONAL, INTENT(in   ) ::   ld4only     ! if .T., do only 4-neighbour comm (ignore corners)
[14021]13      !
[14363]14      INTEGER  ::    ji,  jj,  jk , jl,  jf, jn      ! dummy loop indices
15      INTEGER  ::   ipi, ipj, ipk, ipl, ipf          ! dimension of the input array
[14314]16      INTEGER  ::   ip0i, ip1i, im0i, im1i
17      INTEGER  ::   ip0j, ip1j, im0j, im1j
18      INTEGER  ::   ishti, ishtj, ishti2, ishtj2
[14367]19      INTEGER  ::   iszS, iszR
[14021]20      INTEGER  ::   ierr
[14363]21      INTEGER  ::   ihls, idx
[14314]22      INTEGER  ::   impi_nc
[14343]23      INTEGER  ::   ifill_nfd
[14314]24      INTEGER, DIMENSION(4)  ::   iwewe, issnn
[14363]25      INTEGER, DIMENSION(8)  ::   isizei, ishtSi, ishtRi, ishtPi
26      INTEGER, DIMENSION(8)  ::   isizej, ishtSj, ishtRj, ishtPj
[14314]27      INTEGER, DIMENSION(8)  ::   ifill, iszall
[14834]28      INTEGER, DIMENSION(8)  ::   jnf
[14367]29      INTEGER, DIMENSION(:), ALLOCATABLE  ::   iScnt, iRcnt    ! number of elements to be sent/received
30      INTEGER, DIMENSION(:), ALLOCATABLE  ::   iSdpl, iRdpl    ! displacement in halos arrays
[14314]31      LOGICAL, DIMENSION(8)  ::   llsend, llrecv
32      REAL(PRECISION) ::   zland
[14367]33      LOGICAL  ::   ll4only                                    ! default: 8 neighbourgs
[14021]34      !!----------------------------------------------------------------------
35      !
36      ! ----------------------------------------- !
[14314]37      !     1. local variables initialization     !
[14021]38      ! ----------------------------------------- !
39      !
[14363]40      ipi = SIZE(ptab(1)%pt4d,1)
41      ipj = SIZE(ptab(1)%pt4d,2)
[14349]42      ipk = SIZE(ptab(1)%pt4d,3)
43      ipl = SIZE(ptab(1)%pt4d,4)
44      ipf = kfld
[14021]45      !
46      IF( narea == 1 .AND. numcom == -1 ) CALL mpp_report( cdname, ipk, ipl, ipf, ld_lbc = .TRUE. )
47      !
[14314]48      ! take care of optional parameters
49      !
[14363]50      ihls = nn_hls       ! default definition
51      IF( PRESENT( khls ) )   ihls = khls
52      IF( ihls > n_hlsmax ) THEN
53         WRITE(ctmp1,*) TRIM(cdname), '  is calling lbc_lnk with khls > n_hlsmax : ', khls, '>', n_hlsmax
54         CALL ctl_stop( 'STOP', ctmp1 )
55      ENDIF
56      IF( ipi /= Ni_0+2*ihls ) THEN
57         WRITE(ctmp1,*) TRIM(cdname), ' is calling lbc_lnk with an input array which does not match ihls along i: ', ipi, ihls, Ni_0
58         CALL ctl_stop( 'STOP', ctmp1 )
59      ENDIF
60      IF( ipj /= Nj_0+2*ihls ) THEN
61         WRITE(ctmp1,*) TRIM(cdname), ' is calling lbc_lnk with an input array which does not match ihls along j:', ipj, ihls , Nj_0
62         CALL ctl_stop( 'STOP', ctmp1 )
63      ENDIF
64      !
[14366]65      ll4only = .FALSE.    ! default definition
66      IF( PRESENT(ld4only) )   ll4only = ld4only
[14314]67      !
[14366]68      impi_nc = mpi_nc_com8(ihls)   ! default
69      IF( ll4only )   impi_nc = mpi_nc_com4(ihls)
[14314]70      !
[14021]71      zland = 0._wp                                     ! land filling value: zero by default
72      IF( PRESENT( pfillval ) )   zland = pfillval      ! set land value
73      !
[14314]74      ! define llsend and llrecv: logicals which say if mpi-neibourgs for send or receive exist or not.
75      IF     ( PRESENT(lsend) .AND. PRESENT(lrecv) ) THEN   ! localy defined neighbourgs
[15354]76         CALL ctl_stop( 'STOP', 'mpp_nc_generic+lsend and lrecv not yet implemented')
[14314]77      ELSE IF( PRESENT(lsend) .OR.  PRESENT(lrecv) ) THEN
[14363]78         WRITE(ctmp1,*) TRIM(cdname), '  is calling lbc_lnk with only one of the two arguments lsend or lrecv'
[14314]79         CALL ctl_stop( 'STOP', ctmp1 )
[14366]80      ELSE                                              ! default neighbours
[14363]81         llsend(:) = mpiSnei(ihls,:) >= 0
[14366]82         IF( ll4only )   llsend(5:8) = .FALSE.          ! exclude corners
[14363]83         llrecv(:) = mpiRnei(ihls,:) >= 0
[14366]84         IF( ll4only )   llrecv(5:8) = .FALSE.          ! exclude corners
85      ENDIF
[14021]86      !
[14314]87      ! define ifill: which method should be used to fill each parts (sides+corners) of the halos
88      ! default definition
89      DO jn = 1, 8
90         IF(             llrecv(jn) ) THEN   ;   ifill(jn) = jpfillmpi    ! with an mpi communication
91         ELSEIF(    l_SelfPerio(jn) ) THEN   ;   ifill(jn) = jpfillperio  ! with self-periodicity
92         ELSEIF( PRESENT(kfillmode) ) THEN   ;   ifill(jn) = kfillmode    ! localy defined
93         ELSE                                ;   ifill(jn) = jpfillcst    ! constant value (zland)
[14366]94         ENDIF
[14314]95      END DO
96      ! take care of "indirect self-periodicity" for the corners
97      DO jn = 5, 8
98         IF(.NOT.l_SelfPerio(jn) .AND. l_SelfPerio(jpwe))   ifill(jn) = jpfillnothing   ! no bi-perio but ew-perio: do corners later
99         IF(.NOT.l_SelfPerio(jn) .AND. l_SelfPerio(jpso))   ifill(jn) = jpfillnothing   ! no bi-perio but ns-perio: do corners later
100      END DO
101      ! north fold treatment
[15050]102      IF( l_IdoNFold ) THEN
[14343]103         ifill_nfd = ifill(jpno)             ! if we are here, this means llrecv(jpno) = .false. and l_SelfPerio(jpno) = .false.
104         ifill( (/jpno/) ) = jpfillnothing   ! we do north fold -> do nothing for northern halo
105      ENDIF
[14314]106     
[14343]107      ! We first define the localization and size of the parts of the array that will be sent (s), received (r)
108      ! or used for periodocity (p). The localization is defined as "the bottom left corner - 1" in i and j directions.
109      ! This is a shift that will be applied later in the do loops to pick-up the appropriate part of the array
[14021]110      !
[14363]111      ! all definitions bellow do not refer to N[ij][se]0 so we can use it with any local value of ihls
112      !                   !                       ________________________
113      ip0i =          0   !          im0j = inner |__|________________|__|
114      ip1i =       ihls   !   im1j = inner - halo |  |__|__________|__|  |
115      im1i = ipi-2*ihls   !                       |  |  |          |  |  |
116      im0i = ipi - ihls   !                       |  |  |          |  |  |
117      ip0j =          0   !                       |  |  |          |  |  |
118      ip1j =       ihls   !                       |  |__|__________|__|  |
119      im1j = ipj-2*ihls   !           ip1j = halo |__|__|__________|__|__|
120      im0j = ipj - ihls   !              ip0j = 0 |__|________________|__|
121      !                   !                    ip0i ip1i        im1i im0i
122      !
[14314]123      iwewe(:) = (/ jpwe,jpea,jpwe,jpea /)   ;   issnn(:) = (/ jpso,jpso,jpno,jpno /)
[14363]124      !     sides:     west  east south north      ;   corners: so-we, so-ea, no-we, no-ea
125      isizei(1:4) = (/ ihls, ihls, Ni_0, Ni_0 /)   ;   isizei(5:8) = ihls              ! i- count
126      isizej(1:4) = (/ Nj_0, Nj_0, ihls, ihls /)   ;   isizej(5:8) = ihls              ! j- count
127      ishtSi(1:4) = (/ ip1i, im1i, ip1i, ip1i /)   ;   ishtSi(5:8) = ishtSi( iwewe )   ! i- shift send data
128      ishtSj(1:4) = (/ ip1j, ip1j, ip1j, im1j /)   ;   ishtSj(5:8) = ishtSj( issnn )   ! j- shift send data
129      ishtRi(1:4) = (/ ip0i, im0i, ip1i, ip1i /)   ;   ishtRi(5:8) = ishtRi( iwewe )   ! i- shift received data location
130      ishtRj(1:4) = (/ ip1j, ip1j, ip0j, im0j /)   ;   ishtRj(5:8) = ishtRj( issnn )   ! j- shift received data location
131      ishtPi(1:4) = (/ im1i, ip1i, ip1i, ip1i /)   ;   ishtPi(5:8) = ishtPi( iwewe )   ! i- shift data used for periodicity
132      ishtPj(1:4) = (/ ip1j, ip1j, im1j, ip1j /)   ;   ishtPj(5:8) = ishtPj( issnn )   ! j- shift data used for periodicity
[14021]133      !
[14314]134      ! -------------------------------- !
135      !     2. Prepare MPI exchanges     !
136      ! -------------------------------- !
[14021]137      !
[14314]138      ! Allocate local temporary arrays to be sent/received.
[14367]139      iszS = COUNT( llsend )
140      iszR = COUNT( llrecv )
141      ALLOCATE( iScnt(iszS), iRcnt(iszR), iSdpl(iszS), iRdpl(iszR) )   ! ok if iszS = 0 or iszR = 0
[14314]142      iszall(:) = isizei(:) * isizej(:) * ipk * ipl * ipf
[14367]143      iScnt(:) = PACK( iszall, mask = llsend )                                       ! ok if mask = .false.
144      iRcnt(:) = PACK( iszall, mask = llrecv )
[15301]145      IF( iszS > 0 )   iSdpl(1) = 0
[14367]146      DO jn = 2,iszS
147         iSdpl(jn) = iSdpl(jn-1) + iScnt(jn-1)   ! with _alltoallv: in units of sendtype
[14314]148      END DO
[15301]149      IF( iszR > 0 )   iRdpl(1) = 0
[14367]150      DO jn = 2,iszR
151         iRdpl(jn) = iRdpl(jn-1) + iRcnt(jn-1)   ! with _alltoallv: in units of sendtype
[14314]152      END DO
[14021]153     
[14367]154      ! Allocate buffer arrays to be sent/received if needed
155      iszS = SUM(iszall, mask = llsend)                             ! send buffer size
156      IF( ALLOCATED(BUFFSND) ) THEN
157         IF( SIZE(BUFFSND) < iszS )    DEALLOCATE(BUFFSND)          ! send buffer is too small
158      ENDIF
159      IF( .NOT. ALLOCATED(BUFFSND) )   ALLOCATE( BUFFSND(iszS) )
160      iszR = SUM(iszall, mask = llrecv)                             ! recv buffer size
161      IF( ALLOCATED(BUFFRCV) ) THEN
162         IF( SIZE(BUFFRCV) < iszR )    DEALLOCATE(BUFFRCV)          ! recv buffer is too small
163      ENDIF
164      IF( .NOT. ALLOCATED(BUFFRCV) )   ALLOCATE( BUFFRCV(iszR) )
[14021]165
[14349]166      ! fill sending buffer with ptab(jf)%pt4d
[14021]167      idx = 1
[14314]168      DO jn = 1, 8
169         IF( llsend(jn) ) THEN
[14363]170            ishti = ishtSi(jn)
171            ishtj = ishtSj(jn)
[14314]172            DO jf = 1, ipf  ;  DO jl = 1, ipl  ;  DO jk = 1, ipk  ;  DO jj = 1,isizej(jn)  ;  DO ji = 1,isizei(jn)
[14367]173               BUFFSND(idx) = ptab(jf)%pt4d(ishti+ji,ishtj+jj,jk,jl)
[14021]174               idx = idx + 1
175            END DO   ;   END DO   ;   END DO   ;   END DO   ;   END DO
[14366]176         ENDIF
[14314]177      END DO
[14021]178      !
[14314]179      ! ------------------------------------------------ !
180      !     3. Do all MPI exchanges in 1 unique call     !
181      ! ------------------------------------------------ !
182      !
183      IF( ln_timing ) CALL tic_tac(.TRUE.)
[14367]184      CALL mpi_neighbor_alltoallv (BUFFSND, iScnt, iSdpl, MPI_TYPE, BUFFRCV, iRcnt, iRdpl, MPI_TYPE, impi_nc, ierr)
[14314]185      IF( ln_timing ) CALL tic_tac(.FALSE.)
186      !
187      ! ------------------------- !
188      !     4. Fill all halos     !
189      ! ------------------------- !
190      !
191      idx = 1
[14834]192      ! MPI3 bug fix when domain decomposition has 2 columns/rows
193      IF (jpni .eq. 2) THEN
194         IF (jpnj .eq. 2) THEN
195            jnf(1:8) = (/ 2, 1, 4, 3, 8, 7, 6, 5 /)
196         ELSE
197            jnf(1:8) = (/ 2, 1, 3, 4, 6, 5, 8, 7 /)
198         ENDIF
199      ELSE
200         IF (jpnj .eq. 2) THEN
201            jnf(1:8) = (/ 1, 2, 4, 3, 7, 8, 5, 6 /)
202         ELSE
203            jnf(1:8) = (/ 1, 2, 3, 4, 5, 6, 7, 8 /)
204         ENDIF
205      ENDIF
206
[14314]207      DO jn = 1, 8
[14834]208         ishti = ishtRi(jnf(jn))
209         ishtj = ishtRj(jnf(jn))
210         SELECT CASE ( ifill(jnf(jn)) )
[14021]211         CASE ( jpfillnothing )               ! no filling
[14314]212         CASE ( jpfillmpi   )                 ! fill with data received by MPI
[14834]213            DO jf = 1, ipf  ;  DO jl = 1, ipl  ;  DO jk = 1, ipk  ;  DO jj = 1,isizej(jnf(jn))  ;  DO ji = 1,isizei(jnf(jn))
[14367]214               ptab(jf)%pt4d(ishti+ji,ishtj+jj,jk,jl) = BUFFRCV(idx)
[14021]215               idx = idx + 1
216            END DO   ;   END DO   ;   END DO   ;   END DO   ;   END DO
[14314]217         CASE ( jpfillperio )                 ! use periodicity
[14834]218            ishti2 = ishtPi(jnf(jn))
219            ishtj2 = ishtPj(jnf(jn))
220            DO jf = 1, ipf  ;  DO jl = 1, ipl  ;  DO jk = 1, ipk  ;  DO jj = 1,isizej(jnf(jn))  ;  DO ji = 1,isizei(jnf(jn))
[14349]221               ptab(jf)%pt4d(ishti+ji,ishtj+jj,jk,jl) = ptab(jf)%pt4d(ishti2+ji,ishtj2+jj,jk,jl)
[14021]222            END DO   ;   END DO   ;   END DO   ;   END DO   ;   END DO
223         CASE ( jpfillcopy  )                 ! filling with inner domain values
[14834]224            ishti2 = ishtSi(jnf(jn))
225            ishtj2 = ishtSj(jnf(jn))
226            DO jf = 1, ipf  ;  DO jl = 1, ipl  ;  DO jk = 1, ipk  ;  DO jj = 1,isizej(jnf(jn))  ;  DO ji = 1,isizei(jnf(jn))
[14349]227               ptab(jf)%pt4d(ishti+ji,ishtj+jj,jk,jl) = ptab(jf)%pt4d(ishti2+ji,ishtj2+jj,jk,jl)
[14021]228            END DO   ;   END DO   ;   END DO   ;   END DO   ;   END DO
229         CASE ( jpfillcst   )                 ! filling with constant value
[14834]230            DO jf = 1, ipf  ;  DO jl = 1, ipl  ;  DO jk = 1, ipk  ;  DO jj = 1,isizej(jnf(jn))  ;  DO ji = 1,isizei(jnf(jn))
[14349]231               ptab(jf)%pt4d(ishti+ji,ishtj+jj,jk,jl) = zland
[14021]232            END DO   ;   END DO   ;   END DO   ;   END DO   ;   END DO
233         END SELECT
[14314]234      END DO
[14021]235
[14367]236      DEALLOCATE( iScnt, iRcnt, iSdpl, iRdpl )
237      IF( iszS > jpi*jpj )   DEALLOCATE(BUFFSND)                    ! blocking Send -> can directly deallocate
238      IF( iszR > jpi*jpj )   DEALLOCATE(BUFFRCV)                    ! blocking Recv -> can directly deallocate
[14021]239
[14314]240      ! potential "indirect self-periodicity" for the corners
241      DO jn = 5, 8
242         IF( .NOT. l_SelfPerio(jn) .AND. l_SelfPerio(jpwe)  ) THEN   ! no bi-perio but ew-perio: corners indirect definition
[14363]243            ishti  = ishtRi(jn)
244            ishtj  = ishtRj(jn)
245            ishti2 = ishtPi(jn)   ! use i- shift periodicity
246            ishtj2 = ishtRj(jn)   ! use j- shift recv location: use ew-perio -> ok as filling of the south and north halos now done
[14314]247            DO jf = 1, ipf  ;  DO jl = 1, ipl  ;  DO jk = 1, ipk  ;  DO jj = 1,isizej(jn)  ;  DO ji = 1,isizei(jn)
[14349]248               ptab(jf)%pt4d(ishti+ji,ishtj+jj,jk,jl) = ptab(jf)%pt4d(ishti2+ji,ishtj2+jj,jk,jl)
[14021]249            END DO   ;   END DO   ;   END DO   ;   END DO   ;   END DO
[14314]250         ENDIF
251         IF( .NOT. l_SelfPerio(jn) .AND. l_SelfPerio(jpso)  ) THEN   ! no bi-perio but ns-perio: corners indirect definition
[14363]252            ishti  = ishtRi(jn)
253            ishtj  = ishtRj(jn)
254            ishti2 = ishtRi(jn)   ! use i- shift recv location: use ns-perio -> ok as filling of the west and east halos now done
255            ishtj2 = ishtPj(jn)   ! use j- shift periodicity
[14314]256            DO jf = 1, ipf  ;  DO jl = 1, ipl  ;  DO jk = 1, ipk  ;  DO jj = 1,isizej(jn)  ;  DO ji = 1,isizei(jn)
[14349]257               ptab(jf)%pt4d(ishti+ji,ishtj+jj,jk,jl) = ptab(jf)%pt4d(ishti2+ji,ishtj2+jj,jk,jl)
[14021]258            END DO   ;   END DO   ;   END DO   ;   END DO   ;   END DO
[14314]259         ENDIF
260      END DO
[14021]261      !
262      ! ------------------------------- !
[14314]263      !     5. north fold treatment     !
[14021]264      ! ------------------------------- !
265      !
[15050]266      IF( l_IdoNFold ) THEN
[14363]267         IF( jpni == 1 )  THEN   ;   CALL lbc_nfd( ptab, cd_nat, psgn                  , ihls, ipf )   ! self NFold
268         ELSE                    ;   CALL mpp_nfd( ptab, cd_nat, psgn, ifill_nfd, zland, ihls, ipf )   ! mpi  NFold
[14314]269         ENDIF
[14021]270      ENDIF
[14338]271      !
[14349]272   END SUBROUTINE lbc_lnk_neicoll_/**/PRECISION
[14021]273
Note: See TracBrowser for help on using the repository browser.