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 @ 15104

Last change on this file since 15104 was 15048, checked in by clem, 3 years ago

reduce drastically the number of global communications when using diagnostic outputs. New functions are created: glob_sum_vec and glob_sum_full_vec, which gave a vector as an output. This vector is composed of different summed arrays (such as temperature, salinity etc). Global diagnostics are identical as before. See examples in icedia and diahsb.
  • Property svn:keywords set to Id
File size: 27.3 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)
29   PUBLIC   glob_sum_full ! used in many places (masked with tmask_h, ie only over the halos)
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      !!-----------------------------------------------------------------------
334      !
335      ipi = SIZE(ptab,1)   ! 1st dimension
336      ipj = SIZE(ptab,2)   ! 2nd dimension
337      ipk = SIZE(ptab,3)   ! 3rd dimension
338      !
339      ALLOCATE( ctmp(ipk) )
340      !
341      DO jk = 1, ipk
342         ctmp(jk) = CMPLX( 0.e0, 0.e0, dp )   ! warning ctmp is cumulated
343         DO jj = 1, ipj
344            DO ji = 1, ipi
345               ztmp =  ptab(ji,jj,jk) * tmask_i(ji,jj)
346               CALL DDPDD( CMPLX( ztmp, 0.e0, dp ), ctmp(jk) )
347            END DO
348         END DO
349      END DO
350      CALL mpp_sum( cdname, ctmp(:) )   ! sum over the global domain
351      !
352      ptmp = REAL( ctmp(:), wp )
353      !
354      DEALLOCATE( ctmp )
355      !
356   END FUNCTION glob_sum_vec_3d
357
358   FUNCTION glob_sum_vec_4d( cdname, ptab ) RESULT( ptmp )
359      !!----------------------------------------------------------------------
360      CHARACTER(len=*),  INTENT(in) ::   cdname        ! name of the calling subroutine
361      REAL(wp),          INTENT(in) ::   ptab(:,:,:,:) ! array on which operation is applied
362      REAL(wp), DIMENSION(SIZE(ptab,4)) ::   ptmp
363      !
364      COMPLEX(dp), DIMENSION(:), ALLOCATABLE ::   ctmp
365      REAL(wp)    ::   ztmp
366      INTEGER     ::   ji , jj , jk , jl     ! dummy loop indices
367      INTEGER     ::   ipi, ipj, ipk, ipl    ! dimensions
368      !!-----------------------------------------------------------------------
369      !
370      ipi = SIZE(ptab,1)   ! 1st dimension
371      ipj = SIZE(ptab,2)   ! 2nd dimension
372      ipk = SIZE(ptab,3)   ! 3rd dimension
373      ipl = SIZE(ptab,4)   ! 4th dimension
374      !
375      ALLOCATE( ctmp(ipl) )
376      !
377      DO jl = 1, ipl
378         ctmp(jl) = CMPLX( 0.e0, 0.e0, dp )   ! warning ctmp is cumulated
379         DO jk = 1, ipk
380            DO jj = 1, ipj
381               DO ji = 1, ipi
382                  ztmp =  ptab(ji,jj,jk,jl) * tmask_i(ji,jj)
383                  CALL DDPDD( CMPLX( ztmp, 0.e0, dp ), ctmp(jl) )
384               END DO
385            END DO
386         END DO
387      END DO
388      CALL mpp_sum( cdname, ctmp(:) )   ! sum over the global domain
389      !
390      ptmp = REAL( ctmp(:), wp )
391      !
392      DEALLOCATE( ctmp )
393      !
394   END FUNCTION glob_sum_vec_4d
395   
396   FUNCTION glob_sum_full_vec_3d( cdname, ptab ) RESULT( ptmp )
397      !!----------------------------------------------------------------------
398      CHARACTER(len=*),  INTENT(in) ::   cdname      ! name of the calling subroutine
399      REAL(wp),          INTENT(in) ::   ptab(:,:,:) ! array on which operation is applied
400      REAL(wp), DIMENSION(SIZE(ptab,3)) ::   ptmp
401      !
402      COMPLEX(dp), DIMENSION(:), ALLOCATABLE ::   ctmp
403      REAL(wp)    ::   ztmp
404      INTEGER     ::   ji , jj , jk     ! dummy loop indices
405      INTEGER     ::   ipi, ipj, ipk    ! dimensions
406      !!-----------------------------------------------------------------------
407      !
408      ipi = SIZE(ptab,1)   ! 1st dimension
409      ipj = SIZE(ptab,2)   ! 2nd dimension
410      ipk = SIZE(ptab,3)   ! 3rd dimension
411      !
412      ALLOCATE( ctmp(ipk) )
413      !
414      DO jk = 1, ipk
415         ctmp(jk) = CMPLX( 0.e0, 0.e0, dp )   ! warning ctmp is cumulated
416         DO jj = 1, ipj
417            DO ji = 1, ipi
418               ztmp =  ptab(ji,jj,jk) * tmask_h(ji,jj)
419               CALL DDPDD( CMPLX( ztmp, 0.e0, dp ), ctmp(jk) )
420            END DO
421         END DO
422      END DO
423      CALL mpp_sum( cdname, ctmp(:) )   ! sum over the global domain
424      !
425      ptmp = REAL( ctmp(:), wp )
426      !
427      DEALLOCATE( ctmp )
428      !
429   END FUNCTION glob_sum_full_vec_3d
430
431   FUNCTION glob_sum_full_vec_4d( cdname, ptab ) RESULT( ptmp )
432      !!----------------------------------------------------------------------
433      CHARACTER(len=*),  INTENT(in) ::   cdname        ! name of the calling subroutine
434      REAL(wp),          INTENT(in) ::   ptab(:,:,:,:) ! array on which operation is applied
435      REAL(wp), DIMENSION(SIZE(ptab,4)) ::   ptmp
436      !
437      COMPLEX(dp), DIMENSION(:), ALLOCATABLE ::   ctmp
438      REAL(wp)    ::   ztmp
439      INTEGER     ::   ji , jj , jk , jl     ! dummy loop indices
440      INTEGER     ::   ipi, ipj, ipk, ipl    ! dimensions
441      !!-----------------------------------------------------------------------
442      !
443      ipi = SIZE(ptab,1)   ! 1st dimension
444      ipj = SIZE(ptab,2)   ! 2nd dimension
445      ipk = SIZE(ptab,3)   ! 3rd dimension
446      ipl = SIZE(ptab,4)   ! 4th dimension
447      !
448      ALLOCATE( ctmp(ipl) )
449      !
450      DO jl = 1, ipl
451         ctmp(jl) = CMPLX( 0.e0, 0.e0, dp )   ! warning ctmp is cumulated
452         DO jk = 1, ipk
453            DO jj = 1, ipj
454               DO ji = 1, ipi
455                  ztmp =  ptab(ji,jj,jk,jl) * tmask_h(ji,jj)
456                  CALL DDPDD( CMPLX( ztmp, 0.e0, dp ), ctmp(jl) )
457               END DO
458            END DO
459         END DO
460      END DO
461      CALL mpp_sum( cdname, ctmp(:) )   ! sum over the global domain
462      !
463      ptmp = REAL( ctmp(:), wp )
464      !
465      DEALLOCATE( ctmp )
466      !
467   END FUNCTION glob_sum_full_vec_4d
468
469   SUBROUTINE DDPDD( ydda, yddb )
470      !!----------------------------------------------------------------------
471      !!               ***  ROUTINE DDPDD ***
472      !!
473      !! ** Purpose : Add a scalar element to a sum
474      !!
475      !!
476      !! ** Method  : The code uses the compensated summation with doublet
477      !!              (sum,error) emulated useing complex numbers. ydda is the
478      !!               scalar to add to the summ yddb
479      !!
480      !! ** Action  : This does only work for MPI.
481      !!
482      !! References : Using Acurate Arithmetics to Improve Numerical
483      !!              Reproducibility and Sability in Parallel Applications
484      !!              Yun HE and Chris H. Q. DING, Journal of Supercomputing 18, 259-277, 2001
485      !!----------------------------------------------------------------------
486      COMPLEX(dp), INTENT(in   ) ::   ydda
487      COMPLEX(dp), INTENT(inout) ::   yddb
488      !
489      REAL(dp) :: zerr, zt1, zt2  ! local work variables
490      !!-----------------------------------------------------------------------
491      !
492      ! Compute ydda + yddb using Knuth's trick.
493      zt1  = REAL(ydda) + REAL(yddb)
494      zerr = zt1 - REAL(ydda)
495      zt2  = ( (REAL(yddb) - zerr) + (REAL(ydda) - (zt1 - zerr)) )   &
496         &   + AIMAG(ydda)         + AIMAG(yddb)
497      !
498      ! The result is t1 + t2, after normalization.
499      yddb = CMPLX( zt1 + zt2, zt2 - ((zt1 + zt2) - zt1), wp )
500      !
501   END SUBROUTINE DDPDD
502
503#if defined key_nosignedzero
504   !!----------------------------------------------------------------------
505   !!   'key_nosignedzero'                                         F90 SIGN
506   !!----------------------------------------------------------------------
507
508   FUNCTION SIGN_SCALAR( pa, pb )
509      !!-----------------------------------------------------------------------
510      !!                  ***  FUNCTION SIGN_SCALAR  ***
511      !!
512      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
513      !!-----------------------------------------------------------------------
514      REAL(wp) :: pa,pb          ! input
515      REAL(wp) :: SIGN_SCALAR    ! result
516      !!-----------------------------------------------------------------------
517      IF ( pb >= 0.e0) THEN   ;   SIGN_SCALAR = ABS(pa)
518      ELSE                    ;   SIGN_SCALAR =-ABS(pa)
519      ENDIF
520   END FUNCTION SIGN_SCALAR
521
522
523   FUNCTION SIGN_ARRAY_1D( pa, pb )
524      !!-----------------------------------------------------------------------
525      !!                  ***  FUNCTION SIGN_ARRAY_1D  ***
526      !!
527      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
528      !!-----------------------------------------------------------------------
529      REAL(wp) :: pa,pb(:)                   ! input
530      REAL(wp) :: SIGN_ARRAY_1D(SIZE(pb,1))  ! result
531      !!-----------------------------------------------------------------------
532      WHERE ( pb >= 0.e0 )   ;   SIGN_ARRAY_1D = ABS(pa)
533      ELSEWHERE              ;   SIGN_ARRAY_1D =-ABS(pa)
534      END WHERE
535   END FUNCTION SIGN_ARRAY_1D
536
537
538   FUNCTION SIGN_ARRAY_2D(pa,pb)
539      !!-----------------------------------------------------------------------
540      !!                  ***  FUNCTION SIGN_ARRAY_2D  ***
541      !!
542      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
543      !!-----------------------------------------------------------------------
544      REAL(wp) :: pa,pb(:,:)      ! input
545      REAL(wp) :: SIGN_ARRAY_2D(SIZE(pb,1),SIZE(pb,2))  ! result
546      !!-----------------------------------------------------------------------
547      WHERE ( pb >= 0.e0 )   ;   SIGN_ARRAY_2D = ABS(pa)
548      ELSEWHERE              ;   SIGN_ARRAY_2D =-ABS(pa)
549      END WHERE
550   END FUNCTION SIGN_ARRAY_2D
551
552   FUNCTION SIGN_ARRAY_3D(pa,pb)
553      !!-----------------------------------------------------------------------
554      !!                  ***  FUNCTION SIGN_ARRAY_3D  ***
555      !!
556      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
557      !!-----------------------------------------------------------------------
558      REAL(wp) :: pa,pb(:,:,:)      ! input
559      REAL(wp) :: SIGN_ARRAY_3D(SIZE(pb,1),SIZE(pb,2),SIZE(pb,3))  ! result
560      !!-----------------------------------------------------------------------
561      WHERE ( pb >= 0.e0 )   ;   SIGN_ARRAY_3D = ABS(pa)
562      ELSEWHERE              ;   SIGN_ARRAY_3D =-ABS(pa)
563      END WHERE
564   END FUNCTION SIGN_ARRAY_3D
565
566
567   FUNCTION SIGN_ARRAY_1D_A(pa,pb)
568      !!-----------------------------------------------------------------------
569      !!                  ***  FUNCTION SIGN_ARRAY_1D_A  ***
570      !!
571      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
572      !!-----------------------------------------------------------------------
573      REAL(wp) :: pa(:),pb(:)      ! input
574      REAL(wp) :: SIGN_ARRAY_1D_A(SIZE(pb,1))  ! result
575      !!-----------------------------------------------------------------------
576      WHERE ( pb >= 0.e0 )   ;   SIGN_ARRAY_1D_A = ABS(pa)
577      ELSEWHERE              ;   SIGN_ARRAY_1D_A =-ABS(pa)
578      END WHERE
579   END FUNCTION SIGN_ARRAY_1D_A
580
581
582   FUNCTION SIGN_ARRAY_2D_A(pa,pb)
583      !!-----------------------------------------------------------------------
584      !!                  ***  FUNCTION SIGN_ARRAY_2D_A  ***
585      !!
586      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
587      !!-----------------------------------------------------------------------
588      REAL(wp) :: pa(:,:),pb(:,:)      ! input
589      REAL(wp) :: SIGN_ARRAY_2D_A(SIZE(pb,1),SIZE(pb,2))  ! result
590      !!-----------------------------------------------------------------------
591      WHERE ( pb >= 0.e0 )   ;   SIGN_ARRAY_2D_A = ABS(pa)
592      ELSEWHERE              ;   SIGN_ARRAY_2D_A =-ABS(pa)
593      END WHERE
594   END FUNCTION SIGN_ARRAY_2D_A
595
596
597   FUNCTION SIGN_ARRAY_3D_A(pa,pb)
598      !!-----------------------------------------------------------------------
599      !!                  ***  FUNCTION SIGN_ARRAY_3D_A  ***
600      !!
601      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
602      !!-----------------------------------------------------------------------
603      REAL(wp) :: pa(:,:,:),pb(:,:,:)  ! input
604      REAL(wp) :: SIGN_ARRAY_3D_A(SIZE(pb,1),SIZE(pb,2),SIZE(pb,3)) ! result
605      !!-----------------------------------------------------------------------
606      WHERE ( pb >= 0.e0 )   ;   SIGN_ARRAY_3D_A = ABS(pa)
607      ELSEWHERE              ;   SIGN_ARRAY_3D_A =-ABS(pa)
608      END WHERE
609   END FUNCTION SIGN_ARRAY_3D_A
610
611
612   FUNCTION SIGN_ARRAY_1D_B(pa,pb)
613      !!-----------------------------------------------------------------------
614      !!                  ***  FUNCTION SIGN_ARRAY_1D_B  ***
615      !!
616      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
617      !!-----------------------------------------------------------------------
618      REAL(wp) :: pa(:),pb      ! input
619      REAL(wp) :: SIGN_ARRAY_1D_B(SIZE(pa,1))  ! result
620      !!-----------------------------------------------------------------------
621      IF( pb >= 0.e0 ) THEN   ;   SIGN_ARRAY_1D_B = ABS(pa)
622      ELSE                    ;   SIGN_ARRAY_1D_B =-ABS(pa)
623      ENDIF
624   END FUNCTION SIGN_ARRAY_1D_B
625
626
627   FUNCTION SIGN_ARRAY_2D_B(pa,pb)
628      !!-----------------------------------------------------------------------
629      !!                  ***  FUNCTION SIGN_ARRAY_2D_B  ***
630      !!
631      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
632      !!-----------------------------------------------------------------------
633      REAL(wp) :: pa(:,:),pb      ! input
634      REAL(wp) :: SIGN_ARRAY_2D_B(SIZE(pa,1),SIZE(pa,2))  ! result
635      !!-----------------------------------------------------------------------
636      IF( pb >= 0.e0 ) THEN   ;   SIGN_ARRAY_2D_B = ABS(pa)
637      ELSE                    ;   SIGN_ARRAY_2D_B =-ABS(pa)
638      ENDIF
639   END FUNCTION SIGN_ARRAY_2D_B
640
641
642   FUNCTION SIGN_ARRAY_3D_B(pa,pb)
643      !!-----------------------------------------------------------------------
644      !!                  ***  FUNCTION SIGN_ARRAY_3D_B  ***
645      !!
646      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
647      !!-----------------------------------------------------------------------
648      REAL(wp) :: pa(:,:,:),pb      ! input
649      REAL(wp) :: SIGN_ARRAY_3D_B(SIZE(pa,1),SIZE(pa,2),SIZE(pa,3))  ! result
650      !!-----------------------------------------------------------------------
651      IF( pb >= 0.e0 ) THEN   ;   SIGN_ARRAY_3D_B = ABS(pa)
652      ELSE                    ;   SIGN_ARRAY_3D_B =-ABS(pa)
653      ENDIF
654   END FUNCTION SIGN_ARRAY_3D_B
655#endif
656
657   !!======================================================================
658END MODULE lib_fortran
Note: See TracBrowser for help on using the repository browser.