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

source: NEMO/trunk/src/OCE/LBC/mpp_nfd_generic.h90 @ 13286

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

trunk: merge extra halos branch in trunk, see #2366

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