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/branches/UKMO/dev_r10037_GPU/src/OCE – NEMO

source: NEMO/branches/UKMO/dev_r10037_GPU/src/OCE/lib_fortran.F90 @ 11069

Last change on this file since 11069 was 11069, checked in by andmirek, 5 years ago

ticket #2197 glob sum for parallel architecture

File size: 27.0 KB
RevLine 
[2003]1MODULE lib_fortran
2   !!======================================================================
3   !!                       ***  MODULE  lib_fortran  ***
4   !! Fortran utilities:  includes some low levels fortran functionality
5   !!======================================================================
[2307]6   !! History :  3.2  !  2010-05  (M. Dunphy, R. Benshila)  Original code
[4161]7   !!            3.4  !  2013-06  (C. Rousset)  add glob_min, glob_max
8   !!                                           + 3d dim. of input is fexible (jpk, jpl...)
[7646]9   !!            4.0  !  2016-06  (T. Lovato)  double precision global sum by default
[2003]10   !!----------------------------------------------------------------------
[2307]11
[2003]12   !!----------------------------------------------------------------------
[3764]13   !!   glob_sum    : generic interface for global masked summation over
[2307]14   !!                 the interior domain for 1 or 2 2D or 3D arrays
[3764]15   !!                 it works only for T points
[2307]16   !!   SIGN        : generic interface for SIGN to overwrite f95 behaviour
17   !!                 of intrinsinc sign function
18   !!----------------------------------------------------------------------
[3632]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
[2003]23
24   IMPLICIT NONE
25   PRIVATE
26
[6140]27   PUBLIC   glob_sum      ! used in many places (masked with tmask_i)
[7646]28   PUBLIC   glob_sum_full ! used in many places (masked with tmask_h, ie only over the halos)
[6140]29   PUBLIC   DDPDD         ! also used in closea module
[4161]30   PUBLIC   glob_min, glob_max
[2341]31#if defined key_nosignedzero
[2003]32   PUBLIC SIGN
33#endif
34
35   INTERFACE glob_sum
[3764]36      MODULE PROCEDURE glob_sum_1d, glob_sum_2d, glob_sum_3d, &
37         &             glob_sum_2d_a, glob_sum_3d_a
[2003]38   END INTERFACE
[6140]39   INTERFACE glob_sum_full
40      MODULE PROCEDURE glob_sum_full_2d, glob_sum_full_3d
41   END INTERFACE
[4161]42   INTERFACE glob_min
43      MODULE PROCEDURE glob_min_2d, glob_min_3d,glob_min_2d_a, glob_min_3d_a 
44   END INTERFACE
45   INTERFACE glob_max
46      MODULE PROCEDURE glob_max_2d, glob_max_3d,glob_max_2d_a, glob_max_3d_a 
47   END INTERFACE
[2003]48
[3764]49#if defined key_nosignedzero
[2003]50   INTERFACE SIGN
[2307]51      MODULE PROCEDURE SIGN_SCALAR, SIGN_ARRAY_1D, SIGN_ARRAY_2D, SIGN_ARRAY_3D,   &
[3764]52         &             SIGN_ARRAY_1D_A, SIGN_ARRAY_2D_A, SIGN_ARRAY_3D_A,          &
53         &             SIGN_ARRAY_1D_B, SIGN_ARRAY_2D_B, SIGN_ARRAY_3D_B
[2003]54   END INTERFACE
55#endif
56
[2307]57   !!----------------------------------------------------------------------
[9598]58   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
[3764]59   !! $Id$
[10843]60   !! Software governed by the CeCILL license (see ./LICENSE)
[2307]61   !!----------------------------------------------------------------------
[3764]62CONTAINS
[2003]63
[4161]64   ! --- SUM ---
[3764]65   FUNCTION glob_sum_1d( ptab, kdim )
[6140]66      !!----------------------------------------------------------------------
[3764]67      !!                  ***  FUNCTION  glob_sum_1d ***
68      !!
69      !! ** Purpose : perform a sum in calling DDPDD routine
70      !!----------------------------------------------------------------------
71      INTEGER , INTENT(in) :: kdim
72      REAL(wp), INTENT(in), DIMENSION(kdim) ::   ptab
73      REAL(wp)                              ::   glob_sum_1d   ! global sum
74      !!
75      COMPLEX(wp)::   ctmp
76      REAL(wp)   ::   ztmp
77      INTEGER    ::   ji   ! dummy loop indices
78      !!-----------------------------------------------------------------------
79      !
80      ztmp = 0.e0
81      ctmp = CMPLX( 0.e0, 0.e0, wp )
82      DO ji = 1, kdim
83         ztmp =  ptab(ji)
84         CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
85         END DO
86      IF( lk_mpp )   CALL mpp_sum( ctmp )   ! sum over the global domain
87      glob_sum_1d = REAL(ctmp,wp)
88      !
89   END FUNCTION glob_sum_1d
90
[11069]91   FUNCTION glob_sum_c1d(ptab, kdim)
92      INTEGER, INTENT(IN) :: kdim
93      COMPLEX(KIND = wp), INTENT(IN), DIMENSION(kdim) :: ptab
94      REAL(KIND = wp) :: glob_sum_c1d
95      COMPLEX(KIND = wp) :: ctmp
96      INTEGER :: ji
97
98      ctmp = CMPLX(0.E0, 0.E0, wp)
99
100      DO ji = 1, kdim
101        CALL DDPDD(ptab(ji), ctmp)
102      END DO
103
104      IF (lk_mpp) CALL mpp_sum(ctmp)
105
106      glob_sum_c1d = REAL(ctmp, wp)
107   END FUNCTION glob_sum_c1d
108
[3764]109   FUNCTION glob_sum_2d( ptab )
110      !!----------------------------------------------------------------------
[2307]111      !!                  ***  FUNCTION  glob_sum_2d ***
[2003]112      !!
113      !! ** Purpose : perform a sum in calling DDPDD routine
[2307]114      !!----------------------------------------------------------------------
[4161]115      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptab
116      REAL(wp)                             ::   glob_sum_2d   ! global masked sum
[2003]117      !!
[2307]118      COMPLEX(wp)::   ctmp
119      REAL(wp)   ::   ztmp
120      INTEGER    ::   ji, jj   ! dummy loop indices
[11069]121      COMPLEX(KIND = wp) :: hsum(jpj)
[2307]122      !!-----------------------------------------------------------------------
123      !
124      DO jj = 1, jpj
[11069]125         ctmp = CMPLX( 0.e0, 0.e0, wp )
[2307]126         DO ji =1, jpi
[11069]127            ztmp =  ptab(ji,jj) * tmask_i(ji,jj)
128            CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
[2307]129         END DO
[11069]130         hsum(jj) = ctmp
[2307]131      END DO
[11069]132      glob_sum_2d = glob_sum_c1d(hsum, jpj)
[2307]133      !
[3764]134   END FUNCTION glob_sum_2d
[2307]135
136
[3764]137   FUNCTION glob_sum_3d( ptab )
[2003]138      !!----------------------------------------------------------------------
[2307]139      !!                  ***  FUNCTION  glob_sum_3d ***
140      !!
141      !! ** Purpose : perform a sum on a 3D array in calling DDPDD routine
142      !!----------------------------------------------------------------------
[4161]143      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   ptab
144      REAL(wp)                               ::   glob_sum_3d   ! global masked sum
[2307]145      !!
146      COMPLEX(wp)::   ctmp
147      REAL(wp)   ::   ztmp
148      INTEGER    ::   ji, jj, jk   ! dummy loop indices
[4161]149      INTEGER    ::   ijpk ! local variables: size of ptab
[11069]150      COMPLEX(KIND = wp), allocatable :: hsum(:)
[2307]151      !!-----------------------------------------------------------------------
[2003]152      !
[4161]153      ijpk = SIZE(ptab,3)
[11069]154      ALLOCATE(hsum(ijpk))
[4161]155      !
156      DO jk = 1, ijpk
[11069]157         ctmp = CMPLX( 0.e0, 0.e0, wp )
[2307]158         DO jj = 1, jpj
159            DO ji =1, jpi
160            ztmp =  ptab(ji,jj,jk) * tmask_i(ji,jj)
161            CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
162            END DO
[3764]163         END DO
[11069]164         hsum(jk) = ctmp
[2307]165      END DO
[11069]166      glob_sum_3d = glob_sum_c1d(hsum, ijpk) 
167      DEALLOCATE(hsum)
[2307]168      !
[3764]169   END FUNCTION glob_sum_3d
[2307]170
171
[3764]172   FUNCTION glob_sum_2d_a( ptab1, ptab2 )
[2307]173      !!----------------------------------------------------------------------
174      !!                  ***  FUNCTION  glob_sum_2d_a ***
175      !!
176      !! ** Purpose : perform a sum on two 2D arrays in calling DDPDD routine
177      !!----------------------------------------------------------------------
[4161]178      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptab1, ptab2
179      REAL(wp)                             ::   glob_sum_2d_a   ! global masked sum
[2307]180      !!
181      COMPLEX(wp)::   ctmp
182      REAL(wp)   ::   ztmp
183      INTEGER    ::   ji, jj   ! dummy loop indices
[11069]184      COMPLEX(KIND = wp) :: hsum(jpj)
[2003]185      !!-----------------------------------------------------------------------
[2307]186      !
187      DO jj = 1, jpj
[11069]188         ctmp = CMPLX( 0.e0, 0.e0, wp )
[2003]189         DO ji =1, jpi
[11069]190            ztmp =  ptab1(ji,jj) * tmask_i(ji,jj)
191            CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
192            ztmp =  ptab2(ji,jj) * tmask_i(ji,jj)
193            CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
[2003]194         END DO
[11069]195         hsum(jj) = ctmp
[2003]196      END DO
[11069]197      glob_sum_2d_a = glob_sum_c1d(hsum, jpj)
[2307]198      !
[3764]199   END FUNCTION glob_sum_2d_a
[2003]200
[2307]201
[3764]202   FUNCTION glob_sum_3d_a( ptab1, ptab2 )
[2307]203      !!----------------------------------------------------------------------
204      !!                  ***  FUNCTION  glob_sum_3d_a ***
205      !!
206      !! ** Purpose : perform a sum on two 3D array in calling DDPDD routine
207      !!----------------------------------------------------------------------
[4161]208      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   ptab1, ptab2
209      REAL(wp)                               ::   glob_sum_3d_a   ! global masked sum
[2307]210      !!
211      COMPLEX(wp)::   ctmp
212      REAL(wp)   ::   ztmp
213      INTEGER    ::   ji, jj, jk   ! dummy loop indices
[4161]214      INTEGER    ::   ijpk ! local variables: size of ptab
[11069]215      COMPLEX(KIND = wp), allocatable :: hsum(:)
[2307]216      !!-----------------------------------------------------------------------
217      !
[4161]218      ijpk = SIZE(ptab1,3)
[11069]219      ALLOCATE(hsum(ijpk))
[4161]220      !
221      DO jk = 1, ijpk
[11069]222         ctmp = CMPLX( 0.e0, 0.e0, wp )
[2307]223         DO jj = 1, jpj
[4161]224            DO ji = 1, jpi
225               ztmp =  ptab1(ji,jj,jk) * tmask_i(ji,jj)
226               CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
227               ztmp =  ptab2(ji,jj,jk) * tmask_i(ji,jj)
228               CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
[2307]229            END DO
[4161]230         END DO   
[11069]231         hsum(jk) = ctmp
[2307]232      END DO
[11069]233      glob_sum_3d_a = glob_sum_c1d(hsum, ijpk)
234      DEALLOCATE(hsum)
[2307]235      !
[4161]236   END FUNCTION glob_sum_3d_a   
[2307]237
[6140]238   FUNCTION glob_sum_full_2d( ptab )
239      !!----------------------------------------------------------------------
240      !!                  ***  FUNCTION  glob_sum_full_2d ***
241      !!
242      !! ** Purpose : perform a sum in calling DDPDD routine
243      !!----------------------------------------------------------------------
244      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptab
245      REAL(wp)                             ::   glob_sum_full_2d   ! global sum (nomask)
246      !!
247      COMPLEX(wp)::   ctmp
248      REAL(wp)   ::   ztmp
249      INTEGER    ::   ji, jj   ! dummy loop indices
[11069]250      COMPLEX(KIND = wp) :: hsum(jpj)
[6140]251      !!-----------------------------------------------------------------------
252      !
253      DO jj = 1, jpj
[11069]254         ctmp = CMPLX( 0.e0, 0.e0, wp )
[6140]255         DO ji =1, jpi
[11069]256            ztmp =  ptab(ji,jj) * tmask_h(ji,jj) 
257            CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
[6140]258         END DO
[11069]259         hsum(jj) = ctmp
[6140]260      END DO
[11069]261      glob_sum_full_2d = glob_sum_c1d(hsum, jpj)
[6140]262      !
263   END FUNCTION glob_sum_full_2d
264
265   FUNCTION glob_sum_full_3d( ptab )
266      !!----------------------------------------------------------------------
267      !!                  ***  FUNCTION  glob_sum_full_3d ***
268      !!
269      !! ** Purpose : perform a sum on a 3D array in calling DDPDD routine
270      !!----------------------------------------------------------------------
271      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   ptab
272      REAL(wp)                               ::   glob_sum_full_3d   ! global sum (nomask)
273      !!
274      COMPLEX(wp)::   ctmp
275      REAL(wp)   ::   ztmp
276      INTEGER    ::   ji, jj, jk   ! dummy loop indices
277      INTEGER    ::   ijpk ! local variables: size of ptab
[11069]278      COMPLEX(KIND = wp), allocatable :: hsum(:)
[6140]279      !!-----------------------------------------------------------------------
280      !
281      ijpk = SIZE(ptab,3)
[11069]282      ALLOCATE(hsum(ijpk))
[6140]283      !
284      DO jk = 1, ijpk
[11069]285         ctmp = CMPLX( 0.e0, 0.e0, wp )
[6140]286         DO jj = 1, jpj
287            DO ji =1, jpi
288            ztmp =  ptab(ji,jj,jk) * tmask_h(ji,jj)
289            CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
290            END DO
291         END DO
[11069]292         hsum(jk) = ctmp
[6140]293      END DO
[11069]294      glob_sum_full_3d = glob_sum_c1d(hsum, ijpk)
295      DEALLOCATE(hsum)
[6140]296      !
297   END FUNCTION glob_sum_full_3d
298
[4161]299   ! --- MIN ---
300   FUNCTION glob_min_2d( ptab ) 
301      !!-----------------------------------------------------------------------
302      !!                  ***  FUNCTION  glob_min_2D  ***
303      !!
304      !! ** Purpose : perform a masked min on the inner global domain of a 2D array
305      !!-----------------------------------------------------------------------
306      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptab          ! input 2D array
307      REAL(wp)                             ::   glob_min_2d   ! global masked min
308      !!-----------------------------------------------------------------------
309      !
310      glob_min_2d = MINVAL( ptab(:,:)*tmask_i(:,:) )
311      IF( lk_mpp )   CALL mpp_min( glob_min_2d )
312      !
313   END FUNCTION glob_min_2d
314 
315   FUNCTION glob_min_3d( ptab ) 
316      !!-----------------------------------------------------------------------
317      !!                  ***  FUNCTION  glob_min_3D  ***
318      !!
319      !! ** Purpose : perform a masked min on the inner global domain of a 3D array
320      !!-----------------------------------------------------------------------
321      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   ptab          ! input 3D array
322      REAL(wp)                               ::   glob_min_3d   ! global masked min
323      !!
324      INTEGER :: jk
325      INTEGER :: ijpk ! local variable: size of the 3d dimension of ptab
326      !!-----------------------------------------------------------------------
327      !
328      ijpk = SIZE(ptab,3)
329      !
330      glob_min_3d = MINVAL( ptab(:,:,1)*tmask_i(:,:) )
331      DO jk = 2, ijpk
332         glob_min_3d = MIN( glob_min_3d, MINVAL( ptab(:,:,jk)*tmask_i(:,:) ) )
333      END DO
334      IF( lk_mpp )   CALL mpp_min( glob_min_3d )
335      !
336   END FUNCTION glob_min_3d
337
338
339   FUNCTION glob_min_2d_a( ptab1, ptab2 ) 
340      !!-----------------------------------------------------------------------
341      !!                  ***  FUNCTION  glob_min_2D _a ***
342      !!
343      !! ** Purpose : perform a masked min on the inner global domain of two 2D array
344      !!-----------------------------------------------------------------------
345      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptab1, ptab2    ! input 2D array
346      REAL(wp)            , DIMENSION(2)   ::   glob_min_2d_a   ! global masked min
347      !!-----------------------------------------------------------------------
348      !             
349      glob_min_2d_a(1) = MINVAL( ptab1(:,:)*tmask_i(:,:) )
350      glob_min_2d_a(2) = MINVAL( ptab2(:,:)*tmask_i(:,:) )
351      IF( lk_mpp )   CALL mpp_min( glob_min_2d_a, 2 )
352      !
353   END FUNCTION glob_min_2d_a
354 
355 
356   FUNCTION glob_min_3d_a( ptab1, ptab2 ) 
357      !!-----------------------------------------------------------------------
358      !!                  ***  FUNCTION  glob_min_3D_a ***
359      !!
360      !! ** Purpose : perform a masked min on the inner global domain of two 3D array
361      !!-----------------------------------------------------------------------
362      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   ptab1, ptab2    ! input 3D array
363      REAL(wp)            , DIMENSION(2)     ::   glob_min_3d_a   ! global masked min
364      !!
365      INTEGER :: jk
366      INTEGER :: ijpk ! local variable: size of the 3d dimension of ptab
367      !!-----------------------------------------------------------------------
368      !
369      ijpk = SIZE(ptab1,3)
370      !
371      glob_min_3d_a(1) = MINVAL( ptab1(:,:,1)*tmask_i(:,:) )
372      glob_min_3d_a(2) = MINVAL( ptab2(:,:,1)*tmask_i(:,:) )
373      DO jk = 2, ijpk
374         glob_min_3d_a(1) = MIN( glob_min_3d_a(1), MINVAL( ptab1(:,:,jk)*tmask_i(:,:) ) )
375         glob_min_3d_a(2) = MIN( glob_min_3d_a(2), MINVAL( ptab2(:,:,jk)*tmask_i(:,:) ) )
376      END DO
377      IF( lk_mpp )   CALL mpp_min( glob_min_3d_a, 2 )
378      !
379   END FUNCTION glob_min_3d_a
380
381   ! --- MAX ---
382   FUNCTION glob_max_2d( ptab ) 
383      !!-----------------------------------------------------------------------
384      !!                  ***  FUNCTION  glob_max_2D  ***
385      !!
386      !! ** Purpose : perform a masked max on the inner global domain of a 2D array
387      !!-----------------------------------------------------------------------
388      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptab          ! input 2D array
389      REAL(wp)                             ::   glob_max_2d   ! global masked max
390      !!-----------------------------------------------------------------------
391      !
392      glob_max_2d = MAXVAL( ptab(:,:)*tmask_i(:,:) )
393      IF( lk_mpp )   CALL mpp_max( glob_max_2d )
394      !
395   END FUNCTION glob_max_2d
396 
397   FUNCTION glob_max_3d( ptab ) 
398      !!-----------------------------------------------------------------------
399      !!                  ***  FUNCTION  glob_max_3D  ***
400      !!
401      !! ** Purpose : perform a masked max on the inner global domain of a 3D array
402      !!-----------------------------------------------------------------------
403      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   ptab          ! input 3D array
404      REAL(wp)                               ::   glob_max_3d   ! global masked max
405      !!
406      INTEGER :: jk
407      INTEGER :: ijpk ! local variable: size of the 3d dimension of ptab
408      !!-----------------------------------------------------------------------
409      !
410      ijpk = SIZE(ptab,3)
411      !
412      glob_max_3d = MAXVAL( ptab(:,:,1)*tmask_i(:,:) )
413      DO jk = 2, ijpk
414         glob_max_3d = MAX( glob_max_3d, MAXVAL( ptab(:,:,jk)*tmask_i(:,:) ) )
415      END DO
416      IF( lk_mpp )   CALL mpp_max( glob_max_3d )
417      !
418   END FUNCTION glob_max_3d
419
420
421   FUNCTION glob_max_2d_a( ptab1, ptab2 ) 
422      !!-----------------------------------------------------------------------
423      !!                  ***  FUNCTION  glob_max_2D _a ***
424      !!
425      !! ** Purpose : perform a masked max on the inner global domain of two 2D array
426      !!-----------------------------------------------------------------------
427      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptab1, ptab2    ! input 2D array
428      REAL(wp)            , DIMENSION(2)   ::   glob_max_2d_a   ! global masked max
429      !!-----------------------------------------------------------------------
430      !             
431      glob_max_2d_a(1) = MAXVAL( ptab1(:,:)*tmask_i(:,:) )
432      glob_max_2d_a(2) = MAXVAL( ptab2(:,:)*tmask_i(:,:) )
433      IF( lk_mpp )   CALL mpp_max( glob_max_2d_a, 2 )
434      !
435   END FUNCTION glob_max_2d_a
436 
437 
438   FUNCTION glob_max_3d_a( ptab1, ptab2 ) 
439      !!-----------------------------------------------------------------------
440      !!                  ***  FUNCTION  glob_max_3D_a ***
441      !!
442      !! ** Purpose : perform a masked max on the inner global domain of two 3D array
443      !!-----------------------------------------------------------------------
444      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   ptab1, ptab2    ! input 3D array
445      REAL(wp)            , DIMENSION(2)     ::   glob_max_3d_a   ! global masked max
446      !!
447      INTEGER :: jk
448      INTEGER :: ijpk ! local variable: size of the 3d dimension of ptab
449      !!-----------------------------------------------------------------------
450      !
451      ijpk = SIZE(ptab1,3)
452      !
453      glob_max_3d_a(1) = MAXVAL( ptab1(:,:,1)*tmask_i(:,:) )
454      glob_max_3d_a(2) = MAXVAL( ptab2(:,:,1)*tmask_i(:,:) )
455      DO jk = 2, ijpk
456         glob_max_3d_a(1) = MAX( glob_max_3d_a(1), MAXVAL( ptab1(:,:,jk)*tmask_i(:,:) ) )
457         glob_max_3d_a(2) = MAX( glob_max_3d_a(2), MAXVAL( ptab2(:,:,jk)*tmask_i(:,:) ) )
458      END DO
459      IF( lk_mpp )   CALL mpp_max( glob_max_3d_a, 2 )
460      !
461   END FUNCTION glob_max_3d_a
462
463
[2003]464   SUBROUTINE DDPDD( ydda, yddb )
465      !!----------------------------------------------------------------------
466      !!               ***  ROUTINE DDPDD ***
[3764]467      !!
[2003]468      !! ** Purpose : Add a scalar element to a sum
469      !!
[3764]470      !!
471      !! ** Method  : The code uses the compensated summation with doublet
[2003]472      !!              (sum,error) emulated useing complex numbers. ydda is the
[3764]473      !!               scalar to add to the summ yddb
[2003]474      !!
[3764]475      !! ** Action  : This does only work for MPI.
476      !!
[2003]477      !! References : Using Acurate Arithmetics to Improve Numerical
478      !!              Reproducibility and Sability in Parallel Applications
[3764]479      !!              Yun HE and Chris H. Q. DING, Journal of Supercomputing 18, 259-277, 2001
[2003]480      !!----------------------------------------------------------------------
[2307]481      COMPLEX(wp), INTENT(in   ) ::   ydda
482      COMPLEX(wp), INTENT(inout) ::   yddb
483      !
[2003]484      REAL(wp) :: zerr, zt1, zt2  ! local work variables
[2307]485      !!-----------------------------------------------------------------------
486      !
[2003]487      ! Compute ydda + yddb using Knuth's trick.
[2307]488      zt1  = REAL(ydda) + REAL(yddb)
489      zerr = zt1 - REAL(ydda)
490      zt2  = ( (REAL(yddb) - zerr) + (REAL(ydda) - (zt1 - zerr)) )   &
491         &   + AIMAG(ydda)         + AIMAG(yddb)
492      !
[2003]493      ! The result is t1 + t2, after normalization.
[2307]494      yddb = CMPLX( zt1 + zt2, zt2 - ((zt1 + zt2) - zt1), wp )
495      !
[2003]496   END SUBROUTINE DDPDD
497
498#if defined key_nosignedzero
[2307]499   !!----------------------------------------------------------------------
500   !!   'key_nosignedzero'                                         F90 SIGN
501   !!----------------------------------------------------------------------
[3764]502
[2307]503   FUNCTION SIGN_SCALAR( pa, pb )
[2003]504      !!-----------------------------------------------------------------------
505      !!                  ***  FUNCTION SIGN_SCALAR  ***
506      !!
507      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
508      !!-----------------------------------------------------------------------
509      REAL(wp) :: pa,pb          ! input
[2307]510      REAL(wp) :: SIGN_SCALAR    ! result
511      !!-----------------------------------------------------------------------
512      IF ( pb >= 0.e0) THEN   ;   SIGN_SCALAR = ABS(pa)
513      ELSE                    ;   SIGN_SCALAR =-ABS(pa)
[2003]514      ENDIF
515   END FUNCTION SIGN_SCALAR
516
[2307]517
[3764]518   FUNCTION SIGN_ARRAY_1D( pa, pb )
[2003]519      !!-----------------------------------------------------------------------
520      !!                  ***  FUNCTION SIGN_ARRAY_1D  ***
521      !!
522      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
523      !!-----------------------------------------------------------------------
[2307]524      REAL(wp) :: pa,pb(:)                   ! input
[2003]525      REAL(wp) :: SIGN_ARRAY_1D(SIZE(pb,1))  ! result
[2307]526      !!-----------------------------------------------------------------------
527      WHERE ( pb >= 0.e0 )   ;   SIGN_ARRAY_1D = ABS(pa)
528      ELSEWHERE              ;   SIGN_ARRAY_1D =-ABS(pa)
[2003]529      END WHERE
530   END FUNCTION SIGN_ARRAY_1D
531
[2307]532
[3764]533   FUNCTION SIGN_ARRAY_2D(pa,pb)
[2003]534      !!-----------------------------------------------------------------------
535      !!                  ***  FUNCTION SIGN_ARRAY_2D  ***
536      !!
537      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
538      !!-----------------------------------------------------------------------
539      REAL(wp) :: pa,pb(:,:)      ! input
540      REAL(wp) :: SIGN_ARRAY_2D(SIZE(pb,1),SIZE(pb,2))  ! result
[2307]541      !!-----------------------------------------------------------------------
542      WHERE ( pb >= 0.e0 )   ;   SIGN_ARRAY_2D = ABS(pa)
543      ELSEWHERE              ;   SIGN_ARRAY_2D =-ABS(pa)
[2003]544      END WHERE
545   END FUNCTION SIGN_ARRAY_2D
546
[3764]547   FUNCTION SIGN_ARRAY_3D(pa,pb)
[2003]548      !!-----------------------------------------------------------------------
549      !!                  ***  FUNCTION SIGN_ARRAY_3D  ***
550      !!
551      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
552      !!-----------------------------------------------------------------------
553      REAL(wp) :: pa,pb(:,:,:)      ! input
554      REAL(wp) :: SIGN_ARRAY_3D(SIZE(pb,1),SIZE(pb,2),SIZE(pb,3))  ! result
[2307]555      !!-----------------------------------------------------------------------
556      WHERE ( pb >= 0.e0 )   ;   SIGN_ARRAY_3D = ABS(pa)
557      ELSEWHERE              ;   SIGN_ARRAY_3D =-ABS(pa)
[2003]558      END WHERE
559   END FUNCTION SIGN_ARRAY_3D
560
[2307]561
[3764]562   FUNCTION SIGN_ARRAY_1D_A(pa,pb)
[2003]563      !!-----------------------------------------------------------------------
564      !!                  ***  FUNCTION SIGN_ARRAY_1D_A  ***
565      !!
566      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
567      !!-----------------------------------------------------------------------
568      REAL(wp) :: pa(:),pb(:)      ! input
[2307]569      REAL(wp) :: SIGN_ARRAY_1D_A(SIZE(pb,1))  ! result
570      !!-----------------------------------------------------------------------
571      WHERE ( pb >= 0.e0 )   ;   SIGN_ARRAY_1D_A = ABS(pa)
572      ELSEWHERE              ;   SIGN_ARRAY_1D_A =-ABS(pa)
[2003]573      END WHERE
574   END FUNCTION SIGN_ARRAY_1D_A
575
[2307]576
[3764]577   FUNCTION SIGN_ARRAY_2D_A(pa,pb)
[2003]578      !!-----------------------------------------------------------------------
579      !!                  ***  FUNCTION SIGN_ARRAY_2D_A  ***
580      !!
581      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
582      !!-----------------------------------------------------------------------
583      REAL(wp) :: pa(:,:),pb(:,:)      ! input
584      REAL(wp) :: SIGN_ARRAY_2D_A(SIZE(pb,1),SIZE(pb,2))  ! result
[2307]585      !!-----------------------------------------------------------------------
586      WHERE ( pb >= 0.e0 )   ;   SIGN_ARRAY_2D_A = ABS(pa)
587      ELSEWHERE              ;   SIGN_ARRAY_2D_A =-ABS(pa)
[2003]588      END WHERE
589   END FUNCTION SIGN_ARRAY_2D_A
590
[2307]591
[3764]592   FUNCTION SIGN_ARRAY_3D_A(pa,pb)
[2003]593      !!-----------------------------------------------------------------------
594      !!                  ***  FUNCTION SIGN_ARRAY_3D_A  ***
595      !!
596      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
597      !!-----------------------------------------------------------------------
598      REAL(wp) :: pa(:,:,:),pb(:,:,:)  ! input
599      REAL(wp) :: SIGN_ARRAY_3D_A(SIZE(pb,1),SIZE(pb,2),SIZE(pb,3)) ! result
[2307]600      !!-----------------------------------------------------------------------
601      WHERE ( pb >= 0.e0 )   ;   SIGN_ARRAY_3D_A = ABS(pa)
602      ELSEWHERE              ;   SIGN_ARRAY_3D_A =-ABS(pa)
[2003]603      END WHERE
604   END FUNCTION SIGN_ARRAY_3D_A
605
[2307]606
[3764]607   FUNCTION SIGN_ARRAY_1D_B(pa,pb)
[2003]608      !!-----------------------------------------------------------------------
609      !!                  ***  FUNCTION SIGN_ARRAY_1D_B  ***
610      !!
611      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
612      !!-----------------------------------------------------------------------
613      REAL(wp) :: pa(:),pb      ! input
614      REAL(wp) :: SIGN_ARRAY_1D_B(SIZE(pa,1))  ! result
[2307]615      !!-----------------------------------------------------------------------
616      IF( pb >= 0.e0 ) THEN   ;   SIGN_ARRAY_1D_B = ABS(pa)
617      ELSE                    ;   SIGN_ARRAY_1D_B =-ABS(pa)
[2003]618      ENDIF
619   END FUNCTION SIGN_ARRAY_1D_B
620
[2307]621
[3764]622   FUNCTION SIGN_ARRAY_2D_B(pa,pb)
[2003]623      !!-----------------------------------------------------------------------
624      !!                  ***  FUNCTION SIGN_ARRAY_2D_B  ***
625      !!
626      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
627      !!-----------------------------------------------------------------------
628      REAL(wp) :: pa(:,:),pb      ! input
629      REAL(wp) :: SIGN_ARRAY_2D_B(SIZE(pa,1),SIZE(pa,2))  ! result
[2307]630      !!-----------------------------------------------------------------------
631      IF( pb >= 0.e0 ) THEN   ;   SIGN_ARRAY_2D_B = ABS(pa)
632      ELSE                    ;   SIGN_ARRAY_2D_B =-ABS(pa)
[2003]633      ENDIF
634   END FUNCTION SIGN_ARRAY_2D_B
635
[2307]636
[3764]637   FUNCTION SIGN_ARRAY_3D_B(pa,pb)
[2003]638      !!-----------------------------------------------------------------------
639      !!                  ***  FUNCTION SIGN_ARRAY_3D_B  ***
640      !!
641      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
642      !!-----------------------------------------------------------------------
643      REAL(wp) :: pa(:,:,:),pb      ! input
644      REAL(wp) :: SIGN_ARRAY_3D_B(SIZE(pa,1),SIZE(pa,2),SIZE(pa,3))  ! result
[2307]645      !!-----------------------------------------------------------------------
646      IF( pb >= 0.e0 ) THEN   ;   SIGN_ARRAY_3D_B = ABS(pa)
647      ELSE                    ;   SIGN_ARRAY_3D_B =-ABS(pa)
[2003]648      ENDIF
649   END FUNCTION SIGN_ARRAY_3D_B
650#endif
651
[2307]652   !!======================================================================
[2003]653END MODULE lib_fortran
Note: See TracBrowser for help on using the repository browser.