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

Last change on this file since 13290 was 13290, checked in by smasson, 15 months ago

trunk: mpp_nfd_generic: fix sp/dp compilation issue, following Rachid advise, see #2366

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