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.
mpp_nfd_generic.h90 in NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/OCE/LBC – NEMO

source: NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/OCE/LBC/mpp_nfd_generic.h90 @ 12993

Last change on this file since 12993 was 12993, checked in by smasson, 4 years ago

Extra_Halo: works when removing land subdomain, cleaning/rewriting of mpp_nfd_generic.h90, see #2366

  • Property svn:keywords set to Id
  • Property svn:mime-type set to text/x-fortran
File size: 15.8 KB
Line 
1#if defined MULTI
2#   define NAT_IN(k)                cd_nat(k)   
3#   define SGN_IN(k)                psgn(k)
4#   define F_SIZE(ptab)             kfld
5#   define LBC_ARG                  (jf)
6#   if defined DIM_2d
7#      define ARRAY_TYPE(i,j,k,l,f)    TYPE(PTR_2D)     , INTENT(inout) ::   ptab(f)
8#      define ARRAY_IN(i,j,k,l,f)      ptab(f)%pt2d(i,j)
9#      define K_SIZE(ptab)             1
10#      define L_SIZE(ptab)             1
11#   endif
12#   if defined DIM_3d
13#      define ARRAY_TYPE(i,j,k,l,f)    TYPE(PTR_3D)     , INTENT(inout) ::   ptab(f)
14#      define ARRAY_IN(i,j,k,l,f)      ptab(f)%pt3d(i,j,k)
15#      define K_SIZE(ptab)             SIZE(ptab(1)%pt3d,3)
16#      define L_SIZE(ptab)             1
17#   endif
18#   if defined DIM_4d
19#      define ARRAY_TYPE(i,j,k,l,f)    TYPE(PTR_4D)     , INTENT(inout) ::   ptab(f)
20#      define ARRAY_IN(i,j,k,l,f)      ptab(f)%pt4d(i,j,k,l)
21#      define K_SIZE(ptab)             SIZE(ptab(1)%pt4d,3)
22#      define L_SIZE(ptab)             SIZE(ptab(1)%pt4d,4)
23#   endif
24#else
25!                          !==  IN: ptab is an array  ==!
26#   define ARRAY_TYPE(i,j,k,l,f)    REAL(wp)         , INTENT(inout) ::   ARRAY_IN(i,j,k,l,f)
27#   define NAT_IN(k)                cd_nat
28#   define SGN_IN(k)                psgn
29#   define F_SIZE(ptab)             1
30#   define LBC_ARG
31#   if defined DIM_2d
32#      define ARRAY_IN(i,j,k,l,f)   ptab(i,j)
33#      define K_SIZE(ptab)          1
34#      define L_SIZE(ptab)          1
35#   endif
36#   if defined DIM_3d
37#      define ARRAY_IN(i,j,k,l,f)   ptab(i,j,k)
38#      define K_SIZE(ptab)          SIZE(ptab,3)
39#      define L_SIZE(ptab)          1
40#   endif
41#   if defined DIM_4d
42#      define ARRAY_IN(i,j,k,l,f)   ptab(i,j,k,l)
43#      define K_SIZE(ptab)          SIZE(ptab,3)
44#      define L_SIZE(ptab)          SIZE(ptab,4)
45#   endif
46#endif
47
48   SUBROUTINE ROUTINE_NFD( ptab, cd_nat, psgn, kfillmode, pfillval, kfld )
49      !!----------------------------------------------------------------------
50      ARRAY_TYPE(:,:,:,:,:)   ! array or pointer of arrays on which the boundary condition is applied
51      CHARACTER(len=1) , INTENT(in   ) ::   NAT_IN(:)   ! nature of array grid-points
52      REAL(wp)         , INTENT(in   ) ::   SGN_IN(:)   ! sign used across the north fold boundary
53      INTEGER          , INTENT(in   ) ::   kfillmode   ! filling method for halo over land
54      REAL(wp)         , INTENT(in   ) ::   pfillval    ! background value (used at closed boundaries)
55      INTEGER, OPTIONAL, INTENT(in   ) ::   kfld        ! number of pt3d arrays
56      !
57      LOGICAL  ::   ll_add_line
58      INTEGER  ::   ji,  jj,  jk,  jl, jh, jf, jr   ! dummy loop indices
59      INTEGER  ::   ipi, ipk, ipl, ipf         ! dimension of the input array
60      INTEGER  ::   imigr, iihom, ijhom             ! local integers
61      INTEGER  ::   ierr, ibuffsize, ijpi, iis0, iie0, iilb
62      INTEGER  ::   ijbs, ijbe, ipimax2
63      INTEGER  ::   ij, iproc, ipni
64      INTEGER, DIMENSION (jpmaxngh)       ::   ml_req_nf   ! for mpi_isend when avoiding mpi_allgather
65      INTEGER                             ::   ml_err      ! for mpi_isend when avoiding mpi_allgather
66      INTEGER, DIMENSION(MPI_STATUS_SIZE) ::   ml_stat     ! for mpi_isend when avoiding mpi_allgather
67      !                                                    ! Workspace for message transfers avoiding mpi_allgather
68      INTEGER                             ::   ipj_b       ! sum of lines for all multi fields
69      INTEGER                             ::   ijs, ijb    ! j-counter for send and buffer
70      INTEGER                             ::   i012        ! 0, 1 or 2
71      INTEGER , DIMENSION(:,:)        , ALLOCATABLE ::   jj_s  ! position of sent lines
72      INTEGER , DIMENSION(:,:)        , ALLOCATABLE ::   jj_b  ! position of buffer lines
73      INTEGER , DIMENSION(:)          , ALLOCATABLE ::   ipj_s ! number of sent lines
74      REAL(wp), DIMENSION(:,:,:,:)    , ALLOCATABLE ::   ztabb, ztabr, ztabw  ! buffer, receive and work arrays
75      REAL(wp), DIMENSION(:,:,:,:,:)  , ALLOCATABLE ::   ztab, znorthloc
76      REAL(wp), DIMENSION(:,:,:,:,:,:), ALLOCATABLE ::   znorthgloio
77      !!----------------------------------------------------------------------
78      !
79      ipk = K_SIZE(ptab)   ! 3rd dimension
80      ipl = L_SIZE(ptab)   ! 4th    -
81      ipf = F_SIZE(ptab)   ! 5th    -      use in "multi" case (array of pointers)
82      !
83      IF( l_north_nogather ) THEN      !==  no allgather exchanges  ==!
84
85         !   ---   define number of exchanged lines   ---
86         !
87         ! In theory we should exchange only nn_hls lines.
88         !
89         ! However, some other points are duplicated in the north pole folding:
90         !  - jperio=[34], grid=T : half of the last line (jpiglo/2+2:jpiglo-nn_hls)
91         !  - jperio=[34], grid=U : half of the last line (jpiglo/2+1:jpiglo-nn_hls)
92         !  - jperio=[34], grid=V : all the last line nn_hls+1 and (nn_hls+2:jpiglo-nn_hls)
93         !  - jperio=[34], grid=F : all the last line (nn_hls+1:jpiglo-nn_hls)
94         !  - jperio=[56], grid=T : 2 points of the last line (jpiglo/2+1 and jpglo-nn_hls)
95         !  - jperio=[56], grid=U : no points are duplicated
96         !  - jperio=[56], grid=V : half of the last line (jpiglo/2+1:jpiglo-nn_hls)
97         !  - jperio=[56], grid=F : half of the last line (jpiglo/2+1:jpiglo-nn_hls-1)
98         ! The order of the calculations may differ for these duplicated points (as, for example jj+1 becomes jj-1)
99         ! This explain why these duplicated points may have different values even if they are at the exact same location.
100         ! In consequence, we may want to force the folding on these points by setting l_full_nf_update = .TRUE.
101         ! This is slightly slower but necessary to avoid different values on identical grid points!!
102         !
103         !!!!!!!!!           temporary switch off this optimisation ==> force TRUE           !!!!!!!!
104         !!!!!!!!!  needed to get the same results without agrif and with agrif and no zoom  !!!!!!!!
105         !!!!!!!!!                    I don't know why we must do that...                    !!!!!!!!
106         l_full_nf_update = .TRUE.
107         ! also force it if not restart during the first 2 steps (leap frog?)
108         ll_add_line = l_full_nf_update .OR. ( ncom_stp <= nit000+1 .AND. .NOT. ln_rstart )
109         
110         ALLOCATE(ipj_s(ipf))                ! how many lines do we exchange?
111         IF( ll_add_line ) THEN
112            DO jf = 1, ipf                      ! Loop over the number of arrays to be processed
113               ipj_s(jf) = nn_hls + COUNT( (/ npolj == 3 .OR. npolj == 4 .OR. NAT_IN(jf) == 'V' .OR. NAT_IN(jf) == 'F' /) ) 
114            END DO
115         ELSE
116            ipj_s(:) = nn_hls
117         ENDIF
118         
119         ijpj  = MAXVAL(ipj_s(:))            ! Max 2nd dimension of message transfers (last two j-line only)
120         ipj_b = SUM(   ipj_s(:))            ! Total number of lines to be exchanged
121         ALLOCATE( jj_s(ijpj, ipf), jj_b(ijpj, ipf) )
122
123         ! Index of modifying lines in input
124         ijb = 0
125         DO jf = 1, ipf                      ! Loop over the number of arrays to be processed
126            !
127            SELECT CASE ( npolj )
128            CASE ( 3, 4 )                       ! *  North fold  T-point pivot
129               SELECT CASE ( NAT_IN(jf) )
130               CASE ( 'T', 'W', 'U' )   ;   i012 = 1   ! T-, U-, W-point
131               CASE ( 'V', 'F'      )   ;   i012 = 2   ! V-, F-point
132               END SELECT
133            CASE ( 5, 6 )                       ! *  North fold  F-point pivot
134               SELECT CASE ( NAT_IN(jf) )
135               CASE ( 'T', 'W', 'U' )   ;   i012 = 0   ! T-, U-, W-point
136               CASE ( 'V', 'F'      )   ;   i012 = 1   ! V-, F-point
137               END SELECT
138            END SELECT
139               !
140            DO jj = 1, ipj_s(jf)
141               ijb = ijb + 1
142               jj_b(jj,jf) = ijb
143               jj_s(jj,jf) = jpj - 2*nn_hls + jj - i012
144            END DO
145            !
146         END DO
147         !
148         ALLOCATE( ztabb(jpimax,ipj_b,ipk,ipl) )   ! store all the data to be sent in a buffer array
149         ibuffsize = jpimax * ipj_b * ipk * ipl
150         !
151         DO jf = 1, ipf   ;   DO jl = 1, ipl   ;   DO jk = 1, ipk
152            DO jj = 1, ipj_s(jf)
153               ijb = jj_b(jj,jf)
154               ijs = jj_s(jj,jf)
155               ztabb(    1:jpi   ,ijb,jk,jl) = ARRAY_IN(1:jpi,ijs,jk,jl,jf)
156               ztabb(jpi+1:jpimax,ijb,jk,jl) = 0._wp  ! needed? to avoid sending uninitialized values
157            END DO
158         END DO   ;   END DO   ;   END DO
159         !
160         ! start waiting time measurement
161         IF( ln_timing ) CALL tic_tac(.TRUE.)
162         !
163         ! send the data as soon as possible
164         DO jr = 1, nsndto
165            iproc = nfproc(isendto(jr))
166            IF( iproc /= narea-1 .AND. iproc /= -1 ) THEN
167               CALL mppsend( 5, ztabb, ibuffsize, iproc, ml_req_nf(jr) )
168            ENDIF
169         END DO
170         !
171         ipimax2 = jpimax * jpmaxngh
172         ALLOCATE( ztabw(jpimax,ipj_b,ipk,ipl), ztabr(ipimax2,ipj_b,ipk,ipl) ) 
173         !
174         DO jr = 1, nsndto
175            !
176            ipni  = isendto(jr)
177            iproc = nfproc(ipni)
178            ijpi  = nfjpi (ipni)
179            !
180            IF( ipni ==   1  ) THEN   ;   iis0 =   1             ! domain  left side: as e-w comm already done -> from 1st column
181            ELSE                      ;   iis0 =   1  + nn_hls   ! default: -> from inner domain
182            ENDIF
183            IF( ipni == jpni ) THEN   ;   iie0 = ijpi            ! domain right side: as e-w comm already done -> until last column
184            ELSE                      ;   iie0 = ijpi - nn_hls   ! default: -> until inner domain
185            ENDIF
186            iilb = nfimpp(ipni) - nfimpp(isendto(1))
187            !
188            IF(           iproc == -1 ) THEN   ! No neighbour (land proc that was suppressed)
189               !
190               SELECT CASE ( kfillmode )
191               CASE ( jpfillnothing )               ! no filling
192               CASE ( jpfillcopy    )               ! filling with inner domain values
193                  DO jf = 1, ipf   ;   DO jl = 1, ipl   ;   DO jk = 1, ipk
194                     DO jj = 1, ipj_s(jf)
195                        ijb = jj_b(jj,jf)
196                        ijs = jj_s(jj,jf)
197                        DO ji = iis0, iie0
198                           ztabr(iilb+ji,ijb,jk,jl) = ARRAY_IN(nn_hls+1,ijs,jk,jl,jf)   ! chose to take the 1st iner domain point
199                        END DO
200                     END DO
201                  END DO   ;   END DO   ;   END DO
202               CASE ( jpfillcst     )               ! filling with constant value
203                  DO jl = 1, ipl   ;   DO jk = 1, ipk
204                     DO jj = 1, ipj_b
205                        DO ji = iis0, iie0
206                           ztabr(iilb+ji,jj,jk,jl) = pfillval
207                        END DO
208                     END DO
209                  END DO   ;   END DO
210               END SELECT
211               !
212            ELSE IF( iproc == narea-1 ) THEN   ! get data from myself!
213               !
214               DO jf = 1, ipf   ;   DO jl = 1, ipl  ;   DO jk = 1, ipk
215                  DO jj = 1, ipj_s(jf)
216                     ijb = jj_b(jj,jf)
217                     ijs = jj_s(jj,jf)
218                     DO ji = iis0, iie0
219                        ztabr(iilb+ji,ijb,jk,jl) = ARRAY_IN(ji,ijs,jk,jl,jf)
220                     END DO
221                  END DO
222               END DO   ;   END DO   ;   END DO
223               !
224            ELSE                               ! get data from a neighbour trough communication
225               
226               CALL mpprecv(5, ztabw, ibuffsize, iproc)
227               DO jl = 1, ipl   ;   DO jk = 1, ipk
228                  DO jj = 1, ipj_b
229                     DO ji = iis0, iie0
230                        ztabr(iilb+ji,jj,jk,jl) = ztabw(ji,jj,jk,jl)
231                     END DO
232                  END DO
233               END DO   ;   END DO
234               
235            ENDIF
236         END DO
237         !
238         IF( ln_timing ) CALL tic_tac(.FALSE.)
239         !
240         ! North fold boundary condition
241         !
242         DO jf = 1, ipf
243            ijbs = jj_b(       1 ,jf)
244            ijbe = jj_b(ipj_s(jf),jf)
245            CALL lbc_nfd_nogather( ARRAY_IN(:,:,:,:,jf), ztabr(:,ijbs:ijbe,:,:), cd_nat LBC_ARG, psgn LBC_ARG )
246         END DO
247         !
248         DEALLOCATE( ztabr, ztabw, jj_s, jj_b, ipj_s )
249         !
250         DO jr = 1,nsndto
251            iproc = nfproc(isendto(jr))
252            IF( iproc /= narea-1 .AND. iproc /= -1 ) THEN
253               CALL mpi_wait( ml_req_nf(jr), ml_stat, ml_err )   ! put the wait at the very end just before the deallocate
254            ENDIF
255         END DO
256         DEALLOCATE( ztabb )
257         !
258      ELSE                             !==  allgather exchanges  ==!
259         !
260         ijpj   = 4            ! 2nd dimension of message transfers (last j-lines)
261         !
262         ALLOCATE( znorthloc(jpimax,ijpj,ipk,ipl,ipf) )
263         !
264         DO jf = 1, ipf                ! put in znorthloc the last ijpj j-lines of ptab
265            DO jl = 1, ipl
266               DO jk = 1, ipk
267                  DO jj = jpj - ijpj +1, jpj
268                     ij = jj - jpj + ijpj
269                     znorthloc(1:jpi,ij,jk,jl,jf) = ARRAY_IN(1:jpi,jj,jk,jl,jf)
270                  END DO
271               END DO
272            END DO
273         END DO
274         !
275         ibuffsize = jpimax * ijpj * ipk * ipl * ipf
276         !
277         ALLOCATE( ztab       (jpiglo,ijpj,ipk,ipl,ipf     ) )
278         ALLOCATE( znorthgloio(jpimax,ijpj,ipk,ipl,ipf,jpni) )
279         !
280         ! when some processors of the north fold are suppressed,
281         ! values of ztab* arrays corresponding to these suppressed domain won't be defined
282         ! and we need a default definition to 0.
283         ! a better test should be: a testing if "suppressed land-processors" belongs to the north-pole folding
284         IF ( jpni*jpnj /= jpnij ) ztab(:,:,:,:,:) = 0._wp
285         !
286         ! start waiting time measurement
287         IF( ln_timing ) CALL tic_tac(.TRUE.)
288         CALL MPI_ALLGATHER( znorthloc  , ibuffsize, MPI_DOUBLE_PRECISION,                &
289            &                znorthgloio, ibuffsize, MPI_DOUBLE_PRECISION, ncomm_north, ierr )
290         !
291         ! stop waiting time measurement
292         IF( ln_timing ) CALL tic_tac(.FALSE.)
293         !
294         DO jr = 1, ndim_rank_north         ! recover the global north array
295            iproc = nrank_north(jr) + 1
296            iilb  =  nimppt(iproc)
297            ijpi  =  jpiall(iproc)
298            iis0  = nis0all(iproc)
299            iie0  = nie0all(iproc)
300            IF( iilb            ==      1 )   iis0 = 1      ! e-w boundary already done -> force to take all from 1st column
301            IF( iilb + ijpi - 1 == jpiglo )   iie0 = ijpi   ! e-w boundary already done -> force to take all until last column
302            DO jf = 1, ipf
303               DO jl = 1, ipl
304                  DO jk = 1, ipk
305                     DO jj = 1, ijpj
306                        DO ji = iis0, iie0
307                           ztab(ji+iilb-1,jj,jk,jl,jf) = znorthgloio(ji,jj,jk,jl,jf,jr)
308                        END DO
309                     END DO
310                  END DO
311               END DO
312            END DO
313         END DO
314         DO jf = 1, ipf
315            CALL lbc_nfd( ztab(:,:,:,:,jf), cd_nat LBC_ARG, psgn LBC_ARG )   ! North fold boundary condition
316         END DO
317         !
318         DO jf = 1, ipf
319            DO jl = 1, ipl
320               DO jk = 1, ipk
321                  DO jj = jpj-ijpj+1, jpj             ! Scatter back to ARRAY_IN
322                     ij = jj - jpj + ijpj
323                     DO ji= 1, jpi
324                        ARRAY_IN(ji,jj,jk,jl,jf) = ztab(ji+nimpp-1,ij,jk,jl,jf)
325                     END DO
326                  END DO
327               END DO
328            END DO
329         END DO
330         !
331      !
332         DEALLOCATE( ztab, znorthgloio, znorthloc )
333      ENDIF
334      !
335   END SUBROUTINE ROUTINE_NFD
336
337#undef ARRAY_TYPE
338#undef NAT_IN
339#undef SGN_IN
340#undef ARRAY_IN
341#undef K_SIZE
342#undef L_SIZE
343#undef F_SIZE
344#undef LBC_ARG
Note: See TracBrowser for help on using the repository browser.