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.
lib_fortran.F90 in NEMO/trunk/src/OCE – NEMO

source: NEMO/trunk/src/OCE/lib_fortran.F90 @ 15145

Last change on this file since 15145 was 15145, checked in by smasson, 14 months ago

trunk: pass all sette tests in debug with nn_hls = 2

  • Property svn:keywords set to Id
File size: 29.1 KB
Line 
1MODULE lib_fortran
2   !!======================================================================
3   !!                       ***  MODULE  lib_fortran  ***
4   !! Fortran utilities:  includes some low levels fortran functionality
5   !!======================================================================
6   !! History :  3.2  !  2010-05  (M. Dunphy, R. Benshila)  Original code
7   !!            3.4  !  2013-06  (C. Rousset)  add glob_min, glob_max
8   !!                                           + 3d dim. of input is fexible (jpk, jpl...)
9   !!            4.0  !  2016-06  (T. Lovato)  double precision global sum by default
10   !!----------------------------------------------------------------------
11
12   !!----------------------------------------------------------------------
13   !!   glob_sum    : generic interface for global masked summation over
14   !!                 the interior domain for 1 or 2 2D or 3D arrays
15   !!                 it works only for T points
16   !!   SIGN        : generic interface for SIGN to overwrite f95 behaviour
17   !!                 of intrinsinc sign function
18   !!----------------------------------------------------------------------
19   USE par_oce         ! Ocean parameter
20   USE dom_oce         ! ocean domain
21   USE in_out_manager  ! I/O manager
22   USE lib_mpp         ! distributed memory computing
23   USE lbclnk          ! ocean lateral boundary conditions
24
25   IMPLICIT NONE
26   PRIVATE
27
28   PUBLIC   glob_sum      ! used in many places (masked with tmask_i = ssmask * tmask_h)
29   PUBLIC   glob_sum_full ! used in many places (masked with tmask_h, excluding all duplicated points halos+periodicity)
30   PUBLIC   local_sum     ! used in trcrad, local operation before glob_sum_delay
31   PUBLIC   sum3x3        ! used in trcrad, do a sum over 3x3 boxes
32   PUBLIC   DDPDD         ! also used in closea module
33   PUBLIC   glob_min, glob_max
34   PUBLIC   glob_sum_vec
35   PUBLIC   glob_sum_full_vec
36#if defined key_nosignedzero
37   PUBLIC SIGN
38#endif
39
40   INTERFACE glob_sum
41      MODULE PROCEDURE glob_sum_1d, glob_sum_2d, glob_sum_3d
42   END INTERFACE
43   INTERFACE glob_sum_full
44      MODULE PROCEDURE glob_sum_full_2d, glob_sum_full_3d
45   END INTERFACE
46   INTERFACE local_sum
47      MODULE PROCEDURE local_sum_2d, local_sum_3d
48   END INTERFACE
49   INTERFACE sum3x3
50      MODULE PROCEDURE sum3x3_2d, sum3x3_3d
51   END INTERFACE
52   INTERFACE glob_min
53      MODULE PROCEDURE glob_min_2d, glob_min_3d
54   END INTERFACE
55   INTERFACE glob_max
56      MODULE PROCEDURE glob_max_2d, glob_max_3d
57   END INTERFACE
58   INTERFACE glob_sum_vec
59      MODULE PROCEDURE glob_sum_vec_3d, glob_sum_vec_4d
60   END INTERFACE
61   INTERFACE glob_sum_full_vec
62      MODULE PROCEDURE glob_sum_full_vec_3d, glob_sum_full_vec_4d
63   END INTERFACE
64
65#if defined key_nosignedzero
66   INTERFACE SIGN
67      MODULE PROCEDURE SIGN_SCALAR, SIGN_ARRAY_1D, SIGN_ARRAY_2D, SIGN_ARRAY_3D,   &
68         &             SIGN_ARRAY_1D_A, SIGN_ARRAY_2D_A, SIGN_ARRAY_3D_A,          &
69         &             SIGN_ARRAY_1D_B, SIGN_ARRAY_2D_B, SIGN_ARRAY_3D_B
70   END INTERFACE
71#endif
72
73   !! * Substitutions
74#  include "do_loop_substitute.h90"
75   !!----------------------------------------------------------------------
76   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
77   !! $Id$
78   !! Software governed by the CeCILL license (see ./LICENSE)
79   !!----------------------------------------------------------------------
80CONTAINS
81
82#  define GLOBSUM_CODE
83
84#     define DIM_1d
85#     define FUNCTION_GLOBSUM           glob_sum_1d
86#     include "lib_fortran_generic.h90"
87#     undef FUNCTION_GLOBSUM
88#     undef DIM_1d
89
90#     define DIM_2d
91#     define OPERATION_GLOBSUM
92#     define FUNCTION_GLOBSUM           glob_sum_2d
93#     include "lib_fortran_generic.h90"
94#     undef FUNCTION_GLOBSUM
95#     undef OPERATION_GLOBSUM
96#     define OPERATION_FULL_GLOBSUM
97#     define FUNCTION_GLOBSUM           glob_sum_full_2d
98#     include "lib_fortran_generic.h90"
99#     undef FUNCTION_GLOBSUM
100#     undef OPERATION_FULL_GLOBSUM
101#     undef DIM_2d
102
103#     define DIM_3d
104#     define OPERATION_GLOBSUM
105#     define FUNCTION_GLOBSUM           glob_sum_3d
106#     include "lib_fortran_generic.h90"
107#     undef FUNCTION_GLOBSUM
108#     undef OPERATION_GLOBSUM
109#     define OPERATION_FULL_GLOBSUM
110#     define FUNCTION_GLOBSUM           glob_sum_full_3d
111#     include "lib_fortran_generic.h90"
112#     undef FUNCTION_GLOBSUM
113#     undef OPERATION_FULL_GLOBSUM
114#     undef DIM_3d
115
116#  undef GLOBSUM_CODE
117
118
119#  define GLOBMINMAX_CODE
120
121#     define DIM_2d
122#     define OPERATION_GLOBMIN
123#     define FUNCTION_GLOBMINMAX           glob_min_2d
124#     include "lib_fortran_generic.h90"
125#     undef FUNCTION_GLOBMINMAX
126#     undef OPERATION_GLOBMIN
127#     define OPERATION_GLOBMAX
128#     define FUNCTION_GLOBMINMAX           glob_max_2d
129#     include "lib_fortran_generic.h90"
130#     undef FUNCTION_GLOBMINMAX
131#     undef OPERATION_GLOBMAX
132#     undef DIM_2d
133
134#     define DIM_3d
135#     define OPERATION_GLOBMIN
136#     define FUNCTION_GLOBMINMAX           glob_min_3d
137#     include "lib_fortran_generic.h90"
138#     undef FUNCTION_GLOBMINMAX
139#     undef OPERATION_GLOBMIN
140#     define OPERATION_GLOBMAX
141#     define FUNCTION_GLOBMINMAX           glob_max_3d
142#     include "lib_fortran_generic.h90"
143#     undef FUNCTION_GLOBMINMAX
144#     undef OPERATION_GLOBMAX
145#     undef DIM_3d
146#  undef GLOBMINMAX_CODE
147
148!                          ! FUNCTION local_sum !
149
150   FUNCTION local_sum_2d( ptab )
151      !!----------------------------------------------------------------------
152      REAL(wp),  INTENT(in   ) ::   ptab(:,:) ! array on which operation is applied
153      COMPLEX(dp)              ::  local_sum_2d
154      !
155      !!-----------------------------------------------------------------------
156      !
157      COMPLEX(dp)::   ctmp
158      REAL(wp)   ::   ztmp
159      INTEGER    ::   ji, jj    ! dummy loop indices
160      INTEGER    ::   ipi, ipj  ! dimensions
161      !!-----------------------------------------------------------------------
162      !
163      ipi = SIZE(ptab,1)   ! 1st dimension
164      ipj = SIZE(ptab,2)   ! 2nd dimension
165      !
166      ctmp = CMPLX( 0.e0, 0.e0, wp )   ! warning ctmp is cumulated
167
168      DO jj = 1, ipj
169         DO ji = 1, ipi
170            ztmp =  ptab(ji,jj) * tmask_i(ji,jj)
171            CALL DDPDD( CMPLX( ztmp, 0.e0, dp ), ctmp )
172         END DO
173      END DO
174      !
175      local_sum_2d = ctmp
176       
177   END FUNCTION local_sum_2d
178
179   FUNCTION local_sum_3d( ptab )
180      !!----------------------------------------------------------------------
181      REAL(wp),  INTENT(in   ) ::   ptab(:,:,:) ! array on which operation is applied
182      COMPLEX(dp)              ::  local_sum_3d
183      !
184      !!-----------------------------------------------------------------------
185      !
186      COMPLEX(dp)::   ctmp
187      REAL(wp)   ::   ztmp
188      INTEGER    ::   ji, jj, jk   ! dummy loop indices
189      INTEGER    ::   ipi, ipj, ipk    ! dimensions
190      !!-----------------------------------------------------------------------
191      !
192      ipi = SIZE(ptab,1)   ! 1st dimension
193      ipj = SIZE(ptab,2)   ! 2nd dimension
194      ipk = SIZE(ptab,3)   ! 3rd dimension
195      !
196      ctmp = CMPLX( 0.e0, 0.e0, wp )   ! warning ctmp is cumulated
197
198      DO jk = 1, ipk
199        DO jj = 1, ipj
200          DO ji = 1, ipi
201             ztmp =  ptab(ji,jj,jk) * tmask_i(ji,jj)
202             CALL DDPDD( CMPLX( ztmp, 0.e0, dp ), ctmp )
203          END DO
204        END DO
205      END DO
206      !
207      local_sum_3d = ctmp
208       
209   END FUNCTION local_sum_3d
210
211!                          ! FUNCTION sum3x3 !
212
213   SUBROUTINE sum3x3_2d( p2d )
214      !!-----------------------------------------------------------------------
215      !!                  ***  routine sum3x3_2d  ***
216      !!
217      !! ** Purpose : sum over 3x3 boxes
218      !!----------------------------------------------------------------------
219      REAL(wp), DIMENSION (:,:), INTENT(inout) ::   p2d
220      !
221      INTEGER ::   ji, ji2, jj, jj2     ! dummy loop indices
222      !!----------------------------------------------------------------------
223      !
224      IF( SIZE(p2d,1) /= jpi ) CALL ctl_stop( 'STOP', 'wrong call of sum3x3_2d, the first dimension is not equal to jpi' ) 
225      IF( SIZE(p2d,2) /= jpj ) CALL ctl_stop( 'STOP', 'wrong call of sum3x3_2d, the second dimension is not equal to jpj' ) 
226      !
227      ! work over the whole domain (guarantees all internal cells are set when nn_hls=2)
228      !
229      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
230         IF( MOD(mig(ji), 3) == MOD(nn_hls, 3) .AND.   &              ! 1st bottom left corner always at (Nis0-1, Njs0-1)
231           & MOD(mjg(jj), 3) == MOD(nn_hls, 3)         ) THEN         ! bottom left corner of a 3x3 box
232            ji2 = MIN(mig(ji)+2, jpiglo) - nimpp + 1                  ! right position of the box
233            jj2 = MIN(mjg(jj)+2, jpjglo) - njmpp + 1                  ! upper position of the box
234            IF( ji2 <= jpi .AND. jj2 <= jpj ) THEN                    ! the box is fully included in the local mpi domain
235               p2d(ji:ji2,jj:jj2) = SUM(p2d(ji:ji2,jj:jj2))
236            ENDIF
237         ENDIF
238      END_2D
239      CALL lbc_lnk( 'lib_fortran', p2d, 'T', 1.0_wp )
240      ! no need for 2nd exchange when nn_hls > 1
241      IF( nn_hls == 1 ) THEN
242         IF( mpiRnei(nn_hls,jpwe) > -1 ) THEN   ! 1st column was changed during the previous call to lbc_lnk
243            IF( MOD(mig(    1), 3) == 1 )   &   ! 1st box start at i=1 -> column 1 to 3 correctly computed locally
244               p2d(    1,:) = p2d(    2,:)      ! previous lbc_lnk corrupted column 1 -> put it back using column 2
245            IF( MOD(mig(    1), 3) == 2 )   &   ! 1st box start at i=3 -> column 1 and 2 correctly computed on west neighbourh
246               p2d(    2,:) = p2d(    1,:)      !  previous lbc_lnk fix column 1 -> copy it to column 2
247         ENDIF
248         IF( mpiRnei(nn_hls,jpea) > -1 ) THEN
249            IF( MOD(mig(jpi-2), 3) == 1 )   p2d(  jpi,:) = p2d(jpi-1,:)
250            IF( MOD(mig(jpi-2), 3) == 0 )   p2d(jpi-1,:) = p2d(  jpi,:)
251         ENDIF
252         IF( mpiRnei(nn_hls,jpso) > -1 ) THEN
253            IF( MOD(mjg(    1), 3) == 1 )   p2d(:,    1) = p2d(:,    2)
254            IF( MOD(mjg(    1), 3) == 2 )   p2d(:,    2) = p2d(:,    1)
255         ENDIF
256         IF( mpiRnei(nn_hls,jpno) > -1 ) THEN
257            IF( MOD(mjg(jpj-2), 3) == 1 )   p2d(:,  jpj) = p2d(:,jpj-1)
258            IF( MOD(mjg(jpj-2), 3) == 0 )   p2d(:,jpj-1) = p2d(:,  jpj)
259         ENDIF
260         CALL lbc_lnk( 'lib_fortran', p2d, 'T', 1.0_wp )
261      ENDIF
262
263   END SUBROUTINE sum3x3_2d
264
265   SUBROUTINE sum3x3_3d( p3d )
266      !!-----------------------------------------------------------------------
267      !!                  ***  routine sum3x3_3d  ***
268      !!
269      !! ** Purpose : sum over 3x3 boxes
270      !!----------------------------------------------------------------------
271      REAL(wp), DIMENSION (:,:,:), INTENT(inout) ::   p3d
272      !
273      INTEGER ::   ji, ji2, jj, jj2, jn     ! dummy loop indices
274      INTEGER ::   ipn                      ! Third dimension size
275      !!----------------------------------------------------------------------
276      !
277      IF( SIZE(p3d,1) /= jpi ) CALL ctl_stop( 'STOP', 'wrong call of sum3x3_3d, the first dimension is not equal to jpi' ) 
278      IF( SIZE(p3d,2) /= jpj ) CALL ctl_stop( 'STOP', 'wrong call of sum3x3_3d, the second dimension is not equal to jpj' ) 
279      ipn = SIZE(p3d,3)
280      !
281      DO jn = 1, ipn
282         !
283         ! work over the whole domain (guarantees all internal cells are set when nn_hls=2)
284         !
285         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
286            IF( MOD(mig(ji), 3) == MOD(nn_hls, 3) .AND.   &              ! 1st bottom left corner always at (Nis0-1, Njs0-1)
287              & MOD(mjg(jj), 3) == MOD(nn_hls, 3)         ) THEN         ! bottom left corner of a 3x3 box
288               ji2 = MIN(mig(ji)+2, jpiglo) - nimpp + 1                  ! right position of the box
289               jj2 = MIN(mjg(jj)+2, jpjglo) - njmpp + 1                  ! upper position of the box
290               IF( ji2 <= jpi .AND. jj2 <= jpj ) THEN                    ! the box is fully included in the local mpi domain
291                  p3d(ji:ji2,jj:jj2,jn) = SUM(p3d(ji:ji2,jj:jj2,jn))
292               ENDIF
293            ENDIF
294         END_2D
295      END DO
296      CALL lbc_lnk( 'lib_fortran', p3d, 'T', 1.0_wp )
297      ! no need for 2nd exchange when nn_hls > 1
298      IF( nn_hls == 1 ) THEN
299         IF( mpiRnei(nn_hls,jpwe) > -1 ) THEN    ! 1st column was changed during the previous call to lbc_lnk
300            IF( MOD(mig(    1), 3) == 1 )   &    ! 1st box start at i=1 -> column 1 to 3 correctly computed locally
301               p3d(    1,:,:) = p3d(    2,:,:)   ! previous lbc_lnk corrupted column 1 -> put it back using column 2
302            IF( MOD(mig(    1), 3) == 2 )   &    ! 1st box start at i=3 -> column 1 and 2 correctly computed on west neighbourh
303               p3d(    2,:,:) = p3d(    1,:,:)   !  previous lbc_lnk fix column 1 -> copy it to column 2
304         ENDIF
305         IF( mpiRnei(nn_hls,jpea) > -1 ) THEN
306            IF( MOD(mig(jpi-2), 3) == 1 )   p3d(  jpi,:,:) = p3d(jpi-1,:,:)
307            IF( MOD(mig(jpi-2), 3) == 0 )   p3d(jpi-1,:,:) = p3d(  jpi,:,:)
308         ENDIF
309         IF( mpiRnei(nn_hls,jpso) > -1 ) THEN
310            IF( MOD(mjg(    1), 3) == 1 )   p3d(:,    1,:) = p3d(:,    2,:)
311            IF( MOD(mjg(    1), 3) == 2 )   p3d(:,    2,:) = p3d(:,    1,:)
312         ENDIF
313         IF( mpiRnei(nn_hls,jpno) > -1 ) THEN
314            IF( MOD(mjg(jpj-2), 3) == 1 )   p3d(:,  jpj,:) = p3d(:,jpj-1,:)
315            IF( MOD(mjg(jpj-2), 3) == 0 )   p3d(:,jpj-1,:) = p3d(:,  jpj,:)
316         ENDIF
317         CALL lbc_lnk( 'lib_fortran', p3d, 'T', 1.0_wp )
318      ENDIF
319
320   END SUBROUTINE sum3x3_3d
321
322
323   FUNCTION glob_sum_vec_3d( cdname, ptab ) RESULT( ptmp )
324      !!----------------------------------------------------------------------
325      CHARACTER(len=*),  INTENT(in) ::   cdname      ! name of the calling subroutine
326      REAL(wp),          INTENT(in) ::   ptab(:,:,:) ! array on which operation is applied
327      REAL(wp), DIMENSION(SIZE(ptab,3)) ::   ptmp
328      !
329      COMPLEX(dp), DIMENSION(:), ALLOCATABLE ::   ctmp
330      REAL(wp)    ::   ztmp
331      INTEGER     ::   ji , jj , jk     ! dummy loop indices
332      INTEGER     ::   ipi, ipj, ipk    ! dimensions
333      INTEGER     ::   iis, iie, ijs, ije   ! loop start and end
334      !!-----------------------------------------------------------------------
335      !
336      ipi = SIZE(ptab,1)   ! 1st dimension
337      ipj = SIZE(ptab,2)   ! 2nd dimension
338      ipk = SIZE(ptab,3)   ! 3rd dimension
339      !
340      IF( ipi == jpi .AND. ipj == jpj ) THEN   ! do 2D loop only over the inner domain (-> avoid to use undefined values)
341         iis = Nis0   ;   iie = Nie0
342         ijs = Njs0   ;   ije = Nje0
343      ELSE                                     ! I think we are never in this case...
344         iis = 1   ;   iie = jpi
345         ijs = 1   ;   ije = jpj
346      ENDIF
347      !
348      ALLOCATE( ctmp(ipk) )
349      !
350      DO jk = 1, ipk
351         ctmp(jk) = CMPLX( 0.e0, 0.e0, dp )   ! warning ctmp is cumulated
352         DO jj = ijs, ije
353            DO ji = iis, iie
354               ztmp =  ptab(ji,jj,jk) * tmask_i(ji,jj)
355               CALL DDPDD( CMPLX( ztmp, 0.e0, dp ), ctmp(jk) )
356            END DO
357         END DO
358      END DO
359      CALL mpp_sum( cdname, ctmp(:) )   ! sum over the global domain
360      !
361      ptmp = REAL( ctmp(:), wp )
362      !
363      DEALLOCATE( ctmp )
364      !
365   END FUNCTION glob_sum_vec_3d
366
367   FUNCTION glob_sum_vec_4d( cdname, ptab ) RESULT( ptmp )
368      !!----------------------------------------------------------------------
369      CHARACTER(len=*),  INTENT(in) ::   cdname        ! name of the calling subroutine
370      REAL(wp),          INTENT(in) ::   ptab(:,:,:,:) ! array on which operation is applied
371      REAL(wp), DIMENSION(SIZE(ptab,4)) ::   ptmp
372      !
373      COMPLEX(dp), DIMENSION(:), ALLOCATABLE ::   ctmp
374      REAL(wp)    ::   ztmp
375      INTEGER     ::   ji , jj , jk , jl     ! dummy loop indices
376      INTEGER     ::   ipi, ipj, ipk, ipl    ! dimensions
377      INTEGER     ::   iis, iie, ijs, ije    ! loop start and end
378      !!-----------------------------------------------------------------------
379      !
380      ipi = SIZE(ptab,1)   ! 1st dimension
381      ipj = SIZE(ptab,2)   ! 2nd dimension
382      ipk = SIZE(ptab,3)   ! 3rd dimension
383      ipl = SIZE(ptab,4)   ! 4th dimension
384      !
385      IF( ipi == jpi .AND. ipj == jpj ) THEN   ! do 2D loop only over the inner domain (-> avoid to use undefined values)
386         iis = Nis0   ;   iie = Nie0
387         ijs = Njs0   ;   ije = Nje0
388      ELSE                                     ! I think we are never in this case...
389         iis = 1   ;   iie = jpi
390         ijs = 1   ;   ije = jpj
391      ENDIF
392      !
393      ALLOCATE( ctmp(ipl) )
394      !
395      DO jl = 1, ipl
396         ctmp(jl) = CMPLX( 0.e0, 0.e0, dp )   ! warning ctmp is cumulated
397         DO jk = 1, ipk
398            DO jj = ijs, ije
399               DO ji = iis, iie
400                  ztmp =  ptab(ji,jj,jk,jl) * tmask_i(ji,jj)
401                  CALL DDPDD( CMPLX( ztmp, 0.e0, dp ), ctmp(jl) )
402               END DO
403            END DO
404         END DO
405      END DO
406      CALL mpp_sum( cdname, ctmp(:) )   ! sum over the global domain
407      !
408      ptmp = REAL( ctmp(:), wp )
409      !
410      DEALLOCATE( ctmp )
411      !
412   END FUNCTION glob_sum_vec_4d
413   
414   FUNCTION glob_sum_full_vec_3d( cdname, ptab ) RESULT( ptmp )
415      !!----------------------------------------------------------------------
416      CHARACTER(len=*),  INTENT(in) ::   cdname      ! name of the calling subroutine
417      REAL(wp),          INTENT(in) ::   ptab(:,:,:) ! array on which operation is applied
418      REAL(wp), DIMENSION(SIZE(ptab,3)) ::   ptmp
419      !
420      COMPLEX(dp), DIMENSION(:), ALLOCATABLE ::   ctmp
421      REAL(wp)    ::   ztmp
422      INTEGER     ::   ji , jj , jk     ! dummy loop indices
423      INTEGER     ::   ipi, ipj, ipk    ! dimensions
424      INTEGER     ::   iis, iie, ijs, ije   ! loop start and end
425      !!-----------------------------------------------------------------------
426      !
427      ipi = SIZE(ptab,1)   ! 1st dimension
428      ipj = SIZE(ptab,2)   ! 2nd dimension
429      ipk = SIZE(ptab,3)   ! 3rd dimension
430      !
431      IF( ipi == jpi .AND. ipj == jpj ) THEN   ! do 2D loop only over the inner domain (-> avoid to use undefined values)
432         iis = Nis0   ;   iie = Nie0
433         ijs = Njs0   ;   ije = Nje0
434      ELSE                                     ! I think we are never in this case...
435         iis = 1   ;   iie = jpi
436         ijs = 1   ;   ije = jpj
437      ENDIF
438      !
439      ALLOCATE( ctmp(ipk) )
440      !
441      DO jk = 1, ipk
442         ctmp(jk) = CMPLX( 0.e0, 0.e0, dp )   ! warning ctmp is cumulated
443         DO jj = ijs, ije
444            DO ji = iis, iie
445               ztmp =  ptab(ji,jj,jk) * tmask_h(ji,jj)
446               CALL DDPDD( CMPLX( ztmp, 0.e0, dp ), ctmp(jk) )
447            END DO
448         END DO
449      END DO
450      CALL mpp_sum( cdname, ctmp(:) )   ! sum over the global domain
451      !
452      ptmp = REAL( ctmp(:), wp )
453      !
454      DEALLOCATE( ctmp )
455      !
456   END FUNCTION glob_sum_full_vec_3d
457
458   FUNCTION glob_sum_full_vec_4d( cdname, ptab ) RESULT( ptmp )
459      !!----------------------------------------------------------------------
460      CHARACTER(len=*),  INTENT(in) ::   cdname        ! name of the calling subroutine
461      REAL(wp),          INTENT(in) ::   ptab(:,:,:,:) ! array on which operation is applied
462      REAL(wp), DIMENSION(SIZE(ptab,4)) ::   ptmp
463      !
464      COMPLEX(dp), DIMENSION(:), ALLOCATABLE ::   ctmp
465      REAL(wp)    ::   ztmp
466      INTEGER     ::   ji , jj , jk , jl     ! dummy loop indices
467      INTEGER     ::   ipi, ipj, ipk, ipl    ! dimensions
468      INTEGER     ::   iis, iie, ijs, ije    ! loop start and end
469      !!-----------------------------------------------------------------------
470      !
471      ipi = SIZE(ptab,1)   ! 1st dimension
472      ipj = SIZE(ptab,2)   ! 2nd dimension
473      ipk = SIZE(ptab,3)   ! 3rd dimension
474      ipl = SIZE(ptab,4)   ! 4th dimension
475      !
476      IF( ipi == jpi .AND. ipj == jpj ) THEN   ! do 2D loop only over the inner domain (-> avoid to use undefined values)
477         iis = Nis0   ;   iie = Nie0
478         ijs = Njs0   ;   ije = Nje0
479      ELSE                                     ! I think we are never in this case...
480         iis = 1   ;   iie = jpi
481         ijs = 1   ;   ije = jpj
482      ENDIF
483      !
484      ALLOCATE( ctmp(ipl) )
485      !
486      DO jl = 1, ipl
487         ctmp(jl) = CMPLX( 0.e0, 0.e0, dp )   ! warning ctmp is cumulated
488         DO jk = 1, ipk
489            DO jj = ijs, ije
490               DO ji = iis, iie
491                  ztmp =  ptab(ji,jj,jk,jl) * tmask_h(ji,jj)
492                  CALL DDPDD( CMPLX( ztmp, 0.e0, dp ), ctmp(jl) )
493               END DO
494            END DO
495         END DO
496      END DO
497      CALL mpp_sum( cdname, ctmp(:) )   ! sum over the global domain
498      !
499      ptmp = REAL( ctmp(:), wp )
500      !
501      DEALLOCATE( ctmp )
502      !
503   END FUNCTION glob_sum_full_vec_4d
504
505   SUBROUTINE DDPDD( ydda, yddb )
506      !!----------------------------------------------------------------------
507      !!               ***  ROUTINE DDPDD ***
508      !!
509      !! ** Purpose : Add a scalar element to a sum
510      !!
511      !!
512      !! ** Method  : The code uses the compensated summation with doublet
513      !!              (sum,error) emulated useing complex numbers. ydda is the
514      !!               scalar to add to the summ yddb
515      !!
516      !! ** Action  : This does only work for MPI.
517      !!
518      !! References : Using Acurate Arithmetics to Improve Numerical
519      !!              Reproducibility and Sability in Parallel Applications
520      !!              Yun HE and Chris H. Q. DING, Journal of Supercomputing 18, 259-277, 2001
521      !!----------------------------------------------------------------------
522      COMPLEX(dp), INTENT(in   ) ::   ydda
523      COMPLEX(dp), INTENT(inout) ::   yddb
524      !
525      REAL(dp) :: zerr, zt1, zt2  ! local work variables
526      !!-----------------------------------------------------------------------
527      !
528      ! Compute ydda + yddb using Knuth's trick.
529      zt1  = REAL(ydda) + REAL(yddb)
530      zerr = zt1 - REAL(ydda)
531      zt2  = ( (REAL(yddb) - zerr) + (REAL(ydda) - (zt1 - zerr)) )   &
532         &   + AIMAG(ydda)         + AIMAG(yddb)
533      !
534      ! The result is t1 + t2, after normalization.
535      yddb = CMPLX( zt1 + zt2, zt2 - ((zt1 + zt2) - zt1), wp )
536      !
537   END SUBROUTINE DDPDD
538
539#if defined key_nosignedzero
540   !!----------------------------------------------------------------------
541   !!   'key_nosignedzero'                                         F90 SIGN
542   !!----------------------------------------------------------------------
543
544   FUNCTION SIGN_SCALAR( pa, pb )
545      !!-----------------------------------------------------------------------
546      !!                  ***  FUNCTION SIGN_SCALAR  ***
547      !!
548      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
549      !!-----------------------------------------------------------------------
550      REAL(wp) :: pa,pb          ! input
551      REAL(wp) :: SIGN_SCALAR    ! result
552      !!-----------------------------------------------------------------------
553      IF ( pb >= 0.e0) THEN   ;   SIGN_SCALAR = ABS(pa)
554      ELSE                    ;   SIGN_SCALAR =-ABS(pa)
555      ENDIF
556   END FUNCTION SIGN_SCALAR
557
558
559   FUNCTION SIGN_ARRAY_1D( pa, pb )
560      !!-----------------------------------------------------------------------
561      !!                  ***  FUNCTION SIGN_ARRAY_1D  ***
562      !!
563      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
564      !!-----------------------------------------------------------------------
565      REAL(wp) :: pa,pb(:)                   ! input
566      REAL(wp) :: SIGN_ARRAY_1D(SIZE(pb,1))  ! result
567      !!-----------------------------------------------------------------------
568      WHERE ( pb >= 0.e0 )   ;   SIGN_ARRAY_1D = ABS(pa)
569      ELSEWHERE              ;   SIGN_ARRAY_1D =-ABS(pa)
570      END WHERE
571   END FUNCTION SIGN_ARRAY_1D
572
573
574   FUNCTION SIGN_ARRAY_2D(pa,pb)
575      !!-----------------------------------------------------------------------
576      !!                  ***  FUNCTION SIGN_ARRAY_2D  ***
577      !!
578      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
579      !!-----------------------------------------------------------------------
580      REAL(wp) :: pa,pb(:,:)      ! input
581      REAL(wp) :: SIGN_ARRAY_2D(SIZE(pb,1),SIZE(pb,2))  ! result
582      !!-----------------------------------------------------------------------
583      WHERE ( pb >= 0.e0 )   ;   SIGN_ARRAY_2D = ABS(pa)
584      ELSEWHERE              ;   SIGN_ARRAY_2D =-ABS(pa)
585      END WHERE
586   END FUNCTION SIGN_ARRAY_2D
587
588   FUNCTION SIGN_ARRAY_3D(pa,pb)
589      !!-----------------------------------------------------------------------
590      !!                  ***  FUNCTION SIGN_ARRAY_3D  ***
591      !!
592      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
593      !!-----------------------------------------------------------------------
594      REAL(wp) :: pa,pb(:,:,:)      ! input
595      REAL(wp) :: SIGN_ARRAY_3D(SIZE(pb,1),SIZE(pb,2),SIZE(pb,3))  ! result
596      !!-----------------------------------------------------------------------
597      WHERE ( pb >= 0.e0 )   ;   SIGN_ARRAY_3D = ABS(pa)
598      ELSEWHERE              ;   SIGN_ARRAY_3D =-ABS(pa)
599      END WHERE
600   END FUNCTION SIGN_ARRAY_3D
601
602
603   FUNCTION SIGN_ARRAY_1D_A(pa,pb)
604      !!-----------------------------------------------------------------------
605      !!                  ***  FUNCTION SIGN_ARRAY_1D_A  ***
606      !!
607      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
608      !!-----------------------------------------------------------------------
609      REAL(wp) :: pa(:),pb(:)      ! input
610      REAL(wp) :: SIGN_ARRAY_1D_A(SIZE(pb,1))  ! result
611      !!-----------------------------------------------------------------------
612      WHERE ( pb >= 0.e0 )   ;   SIGN_ARRAY_1D_A = ABS(pa)
613      ELSEWHERE              ;   SIGN_ARRAY_1D_A =-ABS(pa)
614      END WHERE
615   END FUNCTION SIGN_ARRAY_1D_A
616
617
618   FUNCTION SIGN_ARRAY_2D_A(pa,pb)
619      !!-----------------------------------------------------------------------
620      !!                  ***  FUNCTION SIGN_ARRAY_2D_A  ***
621      !!
622      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
623      !!-----------------------------------------------------------------------
624      REAL(wp) :: pa(:,:),pb(:,:)      ! input
625      REAL(wp) :: SIGN_ARRAY_2D_A(SIZE(pb,1),SIZE(pb,2))  ! result
626      !!-----------------------------------------------------------------------
627      WHERE ( pb >= 0.e0 )   ;   SIGN_ARRAY_2D_A = ABS(pa)
628      ELSEWHERE              ;   SIGN_ARRAY_2D_A =-ABS(pa)
629      END WHERE
630   END FUNCTION SIGN_ARRAY_2D_A
631
632
633   FUNCTION SIGN_ARRAY_3D_A(pa,pb)
634      !!-----------------------------------------------------------------------
635      !!                  ***  FUNCTION SIGN_ARRAY_3D_A  ***
636      !!
637      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
638      !!-----------------------------------------------------------------------
639      REAL(wp) :: pa(:,:,:),pb(:,:,:)  ! input
640      REAL(wp) :: SIGN_ARRAY_3D_A(SIZE(pb,1),SIZE(pb,2),SIZE(pb,3)) ! result
641      !!-----------------------------------------------------------------------
642      WHERE ( pb >= 0.e0 )   ;   SIGN_ARRAY_3D_A = ABS(pa)
643      ELSEWHERE              ;   SIGN_ARRAY_3D_A =-ABS(pa)
644      END WHERE
645   END FUNCTION SIGN_ARRAY_3D_A
646
647
648   FUNCTION SIGN_ARRAY_1D_B(pa,pb)
649      !!-----------------------------------------------------------------------
650      !!                  ***  FUNCTION SIGN_ARRAY_1D_B  ***
651      !!
652      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
653      !!-----------------------------------------------------------------------
654      REAL(wp) :: pa(:),pb      ! input
655      REAL(wp) :: SIGN_ARRAY_1D_B(SIZE(pa,1))  ! result
656      !!-----------------------------------------------------------------------
657      IF( pb >= 0.e0 ) THEN   ;   SIGN_ARRAY_1D_B = ABS(pa)
658      ELSE                    ;   SIGN_ARRAY_1D_B =-ABS(pa)
659      ENDIF
660   END FUNCTION SIGN_ARRAY_1D_B
661
662
663   FUNCTION SIGN_ARRAY_2D_B(pa,pb)
664      !!-----------------------------------------------------------------------
665      !!                  ***  FUNCTION SIGN_ARRAY_2D_B  ***
666      !!
667      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
668      !!-----------------------------------------------------------------------
669      REAL(wp) :: pa(:,:),pb      ! input
670      REAL(wp) :: SIGN_ARRAY_2D_B(SIZE(pa,1),SIZE(pa,2))  ! result
671      !!-----------------------------------------------------------------------
672      IF( pb >= 0.e0 ) THEN   ;   SIGN_ARRAY_2D_B = ABS(pa)
673      ELSE                    ;   SIGN_ARRAY_2D_B =-ABS(pa)
674      ENDIF
675   END FUNCTION SIGN_ARRAY_2D_B
676
677
678   FUNCTION SIGN_ARRAY_3D_B(pa,pb)
679      !!-----------------------------------------------------------------------
680      !!                  ***  FUNCTION SIGN_ARRAY_3D_B  ***
681      !!
682      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
683      !!-----------------------------------------------------------------------
684      REAL(wp) :: pa(:,:,:),pb      ! input
685      REAL(wp) :: SIGN_ARRAY_3D_B(SIZE(pa,1),SIZE(pa,2),SIZE(pa,3))  ! result
686      !!-----------------------------------------------------------------------
687      IF( pb >= 0.e0 ) THEN   ;   SIGN_ARRAY_3D_B = ABS(pa)
688      ELSE                    ;   SIGN_ARRAY_3D_B =-ABS(pa)
689      ENDIF
690   END FUNCTION SIGN_ARRAY_3D_B
691#endif
692
693   !!======================================================================
694END MODULE lib_fortran
Note: See TracBrowser for help on using the repository browser.