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 trunk/NEMOGCM/NEMO/OPA_SRC – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/lib_fortran.F90 @ 7881

Last change on this file since 7881 was 7646, checked in by timgraham, 7 years ago

Merge of dev_merge_2016 into trunk. UPDATE TO ARCHFILES NEEDED for XIOS2.
LIM_SRC_s/limrhg.F90 to follow in next commit due to change of kind (I'm unable to do it in this commit).
Merged using the following steps:

1) svn merge --reintegrate svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk .
2) Resolve minor conflicts in sette.sh and namelist_cfg for ORCA2LIM3 (due to a change in trunk after branch was created)
3) svn commit
4) svn switch svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk
5) svn merge svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/2016/dev_merge_2016 .
6) At this stage I checked out a clean copy of the branch to compare against what is about to be committed to the trunk.
6) svn commit #Commit code to the trunk

In this commit I have also reverted a change to Fcheck_archfile.sh which was causing problems on the Paris machine.

  • Property svn:keywords set to Id
File size: 26.4 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   !!----------------------------------------------------------------------
58   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
[3764]59   !! $Id$
[2307]60   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
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
91   FUNCTION glob_sum_2d( ptab )
92      !!----------------------------------------------------------------------
[2307]93      !!                  ***  FUNCTION  glob_sum_2d ***
[2003]94      !!
95      !! ** Purpose : perform a sum in calling DDPDD routine
[2307]96      !!----------------------------------------------------------------------
[4161]97      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptab
98      REAL(wp)                             ::   glob_sum_2d   ! global masked sum
[2003]99      !!
[2307]100      COMPLEX(wp)::   ctmp
101      REAL(wp)   ::   ztmp
102      INTEGER    ::   ji, jj   ! dummy loop indices
103      !!-----------------------------------------------------------------------
104      !
105      ztmp = 0.e0
106      ctmp = CMPLX( 0.e0, 0.e0, wp )
107      DO jj = 1, jpj
108         DO ji =1, jpi
109         ztmp =  ptab(ji,jj) * tmask_i(ji,jj)
110         CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
111         END DO
112      END DO
113      IF( lk_mpp )   CALL mpp_sum( ctmp )   ! sum over the global domain
[3294]114      glob_sum_2d = REAL(ctmp,wp)
[2307]115      !
[3764]116   END FUNCTION glob_sum_2d
[2307]117
118
[3764]119   FUNCTION glob_sum_3d( ptab )
[2003]120      !!----------------------------------------------------------------------
[2307]121      !!                  ***  FUNCTION  glob_sum_3d ***
122      !!
123      !! ** Purpose : perform a sum on a 3D array in calling DDPDD routine
124      !!----------------------------------------------------------------------
[4161]125      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   ptab
126      REAL(wp)                               ::   glob_sum_3d   ! global masked sum
[2307]127      !!
128      COMPLEX(wp)::   ctmp
129      REAL(wp)   ::   ztmp
130      INTEGER    ::   ji, jj, jk   ! dummy loop indices
[4161]131      INTEGER    ::   ijpk ! local variables: size of ptab
[2307]132      !!-----------------------------------------------------------------------
[2003]133      !
[4161]134      ijpk = SIZE(ptab,3)
135      !
[2307]136      ztmp = 0.e0
137      ctmp = CMPLX( 0.e0, 0.e0, wp )
[4161]138      DO jk = 1, ijpk
[2307]139         DO jj = 1, jpj
140            DO ji =1, jpi
141            ztmp =  ptab(ji,jj,jk) * tmask_i(ji,jj)
142            CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
143            END DO
[3764]144         END DO
[2307]145      END DO
146      IF( lk_mpp )   CALL mpp_sum( ctmp )   ! sum over the global domain
[3294]147      glob_sum_3d = REAL(ctmp,wp)
[2307]148      !
[3764]149   END FUNCTION glob_sum_3d
[2307]150
151
[3764]152   FUNCTION glob_sum_2d_a( ptab1, ptab2 )
[2307]153      !!----------------------------------------------------------------------
154      !!                  ***  FUNCTION  glob_sum_2d_a ***
155      !!
156      !! ** Purpose : perform a sum on two 2D arrays in calling DDPDD routine
157      !!----------------------------------------------------------------------
[4161]158      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptab1, ptab2
159      REAL(wp)                             ::   glob_sum_2d_a   ! global masked sum
[2307]160      !!
161      COMPLEX(wp)::   ctmp
162      REAL(wp)   ::   ztmp
163      INTEGER    ::   ji, jj   ! dummy loop indices
[2003]164      !!-----------------------------------------------------------------------
[2307]165      !
[2003]166      ztmp = 0.e0
[2307]167      ctmp = CMPLX( 0.e0, 0.e0, wp )
168      DO jj = 1, jpj
[2003]169         DO ji =1, jpi
[2307]170         ztmp =  ptab1(ji,jj) * tmask_i(ji,jj)
171         CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
172         ztmp =  ptab2(ji,jj) * tmask_i(ji,jj)
173         CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
[2003]174         END DO
175      END DO
176      IF( lk_mpp )   CALL mpp_sum( ctmp )   ! sum over the global domain
[3294]177      glob_sum_2d_a = REAL(ctmp,wp)
[2307]178      !
[3764]179   END FUNCTION glob_sum_2d_a
[2003]180
[2307]181
[3764]182   FUNCTION glob_sum_3d_a( ptab1, ptab2 )
[2307]183      !!----------------------------------------------------------------------
184      !!                  ***  FUNCTION  glob_sum_3d_a ***
185      !!
186      !! ** Purpose : perform a sum on two 3D array in calling DDPDD routine
187      !!----------------------------------------------------------------------
[4161]188      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   ptab1, ptab2
189      REAL(wp)                               ::   glob_sum_3d_a   ! global masked sum
[2307]190      !!
191      COMPLEX(wp)::   ctmp
192      REAL(wp)   ::   ztmp
193      INTEGER    ::   ji, jj, jk   ! dummy loop indices
[4161]194      INTEGER    ::   ijpk ! local variables: size of ptab
[2307]195      !!-----------------------------------------------------------------------
196      !
[4161]197      ijpk = SIZE(ptab1,3)
198      !
[2307]199      ztmp = 0.e0
200      ctmp = CMPLX( 0.e0, 0.e0, wp )
[4161]201      DO jk = 1, ijpk
[2307]202         DO jj = 1, jpj
[4161]203            DO ji = 1, jpi
204               ztmp =  ptab1(ji,jj,jk) * tmask_i(ji,jj)
205               CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
206               ztmp =  ptab2(ji,jj,jk) * tmask_i(ji,jj)
207               CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
[2307]208            END DO
[4161]209         END DO   
[2307]210      END DO
211      IF( lk_mpp )   CALL mpp_sum( ctmp )   ! sum over the global domain
[3294]212      glob_sum_3d_a = REAL(ctmp,wp)
[2307]213      !
[4161]214   END FUNCTION glob_sum_3d_a   
[2307]215
[6140]216   FUNCTION glob_sum_full_2d( ptab )
217      !!----------------------------------------------------------------------
218      !!                  ***  FUNCTION  glob_sum_full_2d ***
219      !!
220      !! ** Purpose : perform a sum in calling DDPDD routine
221      !!----------------------------------------------------------------------
222      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptab
223      REAL(wp)                             ::   glob_sum_full_2d   ! global sum (nomask)
224      !!
225      COMPLEX(wp)::   ctmp
226      REAL(wp)   ::   ztmp
227      INTEGER    ::   ji, jj   ! dummy loop indices
228      !!-----------------------------------------------------------------------
229      !
230      ztmp = 0.e0
231      ctmp = CMPLX( 0.e0, 0.e0, wp )
232      DO jj = 1, jpj
233         DO ji =1, jpi
234         ztmp =  ptab(ji,jj) * tmask_h(ji,jj) 
235         CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
236         END DO
237      END DO
238      IF( lk_mpp )   CALL mpp_sum( ctmp )   ! sum over the global domain
239      glob_sum_full_2d = REAL(ctmp,wp)
240      !
241   END FUNCTION glob_sum_full_2d
242
243   FUNCTION glob_sum_full_3d( ptab )
244      !!----------------------------------------------------------------------
245      !!                  ***  FUNCTION  glob_sum_full_3d ***
246      !!
247      !! ** Purpose : perform a sum on a 3D array in calling DDPDD routine
248      !!----------------------------------------------------------------------
249      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   ptab
250      REAL(wp)                               ::   glob_sum_full_3d   ! global sum (nomask)
251      !!
252      COMPLEX(wp)::   ctmp
253      REAL(wp)   ::   ztmp
254      INTEGER    ::   ji, jj, jk   ! dummy loop indices
255      INTEGER    ::   ijpk ! local variables: size of ptab
256      !!-----------------------------------------------------------------------
257      !
258      ijpk = SIZE(ptab,3)
259      !
260      ztmp = 0.e0
261      ctmp = CMPLX( 0.e0, 0.e0, wp )
262      DO jk = 1, ijpk
263         DO jj = 1, jpj
264            DO ji =1, jpi
265            ztmp =  ptab(ji,jj,jk) * tmask_h(ji,jj)
266            CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
267            END DO
268         END DO
269      END DO
270      IF( lk_mpp )   CALL mpp_sum( ctmp )   ! sum over the global domain
271      glob_sum_full_3d = REAL(ctmp,wp)
272      !
273   END FUNCTION glob_sum_full_3d
274
[4161]275   ! --- MIN ---
276   FUNCTION glob_min_2d( ptab ) 
277      !!-----------------------------------------------------------------------
278      !!                  ***  FUNCTION  glob_min_2D  ***
279      !!
280      !! ** Purpose : perform a masked min on the inner global domain of a 2D array
281      !!-----------------------------------------------------------------------
282      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptab          ! input 2D array
283      REAL(wp)                             ::   glob_min_2d   ! global masked min
284      !!-----------------------------------------------------------------------
285      !
286      glob_min_2d = MINVAL( ptab(:,:)*tmask_i(:,:) )
287      IF( lk_mpp )   CALL mpp_min( glob_min_2d )
288      !
289   END FUNCTION glob_min_2d
290 
291   FUNCTION glob_min_3d( ptab ) 
292      !!-----------------------------------------------------------------------
293      !!                  ***  FUNCTION  glob_min_3D  ***
294      !!
295      !! ** Purpose : perform a masked min on the inner global domain of a 3D array
296      !!-----------------------------------------------------------------------
297      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   ptab          ! input 3D array
298      REAL(wp)                               ::   glob_min_3d   ! global masked min
299      !!
300      INTEGER :: jk
301      INTEGER :: ijpk ! local variable: size of the 3d dimension of ptab
302      !!-----------------------------------------------------------------------
303      !
304      ijpk = SIZE(ptab,3)
305      !
306      glob_min_3d = MINVAL( ptab(:,:,1)*tmask_i(:,:) )
307      DO jk = 2, ijpk
308         glob_min_3d = MIN( glob_min_3d, MINVAL( ptab(:,:,jk)*tmask_i(:,:) ) )
309      END DO
310      IF( lk_mpp )   CALL mpp_min( glob_min_3d )
311      !
312   END FUNCTION glob_min_3d
313
314
315   FUNCTION glob_min_2d_a( ptab1, ptab2 ) 
316      !!-----------------------------------------------------------------------
317      !!                  ***  FUNCTION  glob_min_2D _a ***
318      !!
319      !! ** Purpose : perform a masked min on the inner global domain of two 2D array
320      !!-----------------------------------------------------------------------
321      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptab1, ptab2    ! input 2D array
322      REAL(wp)            , DIMENSION(2)   ::   glob_min_2d_a   ! global masked min
323      !!-----------------------------------------------------------------------
324      !             
325      glob_min_2d_a(1) = MINVAL( ptab1(:,:)*tmask_i(:,:) )
326      glob_min_2d_a(2) = MINVAL( ptab2(:,:)*tmask_i(:,:) )
327      IF( lk_mpp )   CALL mpp_min( glob_min_2d_a, 2 )
328      !
329   END FUNCTION glob_min_2d_a
330 
331 
332   FUNCTION glob_min_3d_a( ptab1, ptab2 ) 
333      !!-----------------------------------------------------------------------
334      !!                  ***  FUNCTION  glob_min_3D_a ***
335      !!
336      !! ** Purpose : perform a masked min on the inner global domain of two 3D array
337      !!-----------------------------------------------------------------------
338      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   ptab1, ptab2    ! input 3D array
339      REAL(wp)            , DIMENSION(2)     ::   glob_min_3d_a   ! global masked min
340      !!
341      INTEGER :: jk
342      INTEGER :: ijpk ! local variable: size of the 3d dimension of ptab
343      !!-----------------------------------------------------------------------
344      !
345      ijpk = SIZE(ptab1,3)
346      !
347      glob_min_3d_a(1) = MINVAL( ptab1(:,:,1)*tmask_i(:,:) )
348      glob_min_3d_a(2) = MINVAL( ptab2(:,:,1)*tmask_i(:,:) )
349      DO jk = 2, ijpk
350         glob_min_3d_a(1) = MIN( glob_min_3d_a(1), MINVAL( ptab1(:,:,jk)*tmask_i(:,:) ) )
351         glob_min_3d_a(2) = MIN( glob_min_3d_a(2), MINVAL( ptab2(:,:,jk)*tmask_i(:,:) ) )
352      END DO
353      IF( lk_mpp )   CALL mpp_min( glob_min_3d_a, 2 )
354      !
355   END FUNCTION glob_min_3d_a
356
357   ! --- MAX ---
358   FUNCTION glob_max_2d( ptab ) 
359      !!-----------------------------------------------------------------------
360      !!                  ***  FUNCTION  glob_max_2D  ***
361      !!
362      !! ** Purpose : perform a masked max on the inner global domain of a 2D array
363      !!-----------------------------------------------------------------------
364      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptab          ! input 2D array
365      REAL(wp)                             ::   glob_max_2d   ! global masked max
366      !!-----------------------------------------------------------------------
367      !
368      glob_max_2d = MAXVAL( ptab(:,:)*tmask_i(:,:) )
369      IF( lk_mpp )   CALL mpp_max( glob_max_2d )
370      !
371   END FUNCTION glob_max_2d
372 
373   FUNCTION glob_max_3d( ptab ) 
374      !!-----------------------------------------------------------------------
375      !!                  ***  FUNCTION  glob_max_3D  ***
376      !!
377      !! ** Purpose : perform a masked max on the inner global domain of a 3D array
378      !!-----------------------------------------------------------------------
379      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   ptab          ! input 3D array
380      REAL(wp)                               ::   glob_max_3d   ! global masked max
381      !!
382      INTEGER :: jk
383      INTEGER :: ijpk ! local variable: size of the 3d dimension of ptab
384      !!-----------------------------------------------------------------------
385      !
386      ijpk = SIZE(ptab,3)
387      !
388      glob_max_3d = MAXVAL( ptab(:,:,1)*tmask_i(:,:) )
389      DO jk = 2, ijpk
390         glob_max_3d = MAX( glob_max_3d, MAXVAL( ptab(:,:,jk)*tmask_i(:,:) ) )
391      END DO
392      IF( lk_mpp )   CALL mpp_max( glob_max_3d )
393      !
394   END FUNCTION glob_max_3d
395
396
397   FUNCTION glob_max_2d_a( ptab1, ptab2 ) 
398      !!-----------------------------------------------------------------------
399      !!                  ***  FUNCTION  glob_max_2D _a ***
400      !!
401      !! ** Purpose : perform a masked max on the inner global domain of two 2D array
402      !!-----------------------------------------------------------------------
403      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptab1, ptab2    ! input 2D array
404      REAL(wp)            , DIMENSION(2)   ::   glob_max_2d_a   ! global masked max
405      !!-----------------------------------------------------------------------
406      !             
407      glob_max_2d_a(1) = MAXVAL( ptab1(:,:)*tmask_i(:,:) )
408      glob_max_2d_a(2) = MAXVAL( ptab2(:,:)*tmask_i(:,:) )
409      IF( lk_mpp )   CALL mpp_max( glob_max_2d_a, 2 )
410      !
411   END FUNCTION glob_max_2d_a
412 
413 
414   FUNCTION glob_max_3d_a( ptab1, ptab2 ) 
415      !!-----------------------------------------------------------------------
416      !!                  ***  FUNCTION  glob_max_3D_a ***
417      !!
418      !! ** Purpose : perform a masked max on the inner global domain of two 3D array
419      !!-----------------------------------------------------------------------
420      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   ptab1, ptab2    ! input 3D array
421      REAL(wp)            , DIMENSION(2)     ::   glob_max_3d_a   ! global masked max
422      !!
423      INTEGER :: jk
424      INTEGER :: ijpk ! local variable: size of the 3d dimension of ptab
425      !!-----------------------------------------------------------------------
426      !
427      ijpk = SIZE(ptab1,3)
428      !
429      glob_max_3d_a(1) = MAXVAL( ptab1(:,:,1)*tmask_i(:,:) )
430      glob_max_3d_a(2) = MAXVAL( ptab2(:,:,1)*tmask_i(:,:) )
431      DO jk = 2, ijpk
432         glob_max_3d_a(1) = MAX( glob_max_3d_a(1), MAXVAL( ptab1(:,:,jk)*tmask_i(:,:) ) )
433         glob_max_3d_a(2) = MAX( glob_max_3d_a(2), MAXVAL( ptab2(:,:,jk)*tmask_i(:,:) ) )
434      END DO
435      IF( lk_mpp )   CALL mpp_max( glob_max_3d_a, 2 )
436      !
437   END FUNCTION glob_max_3d_a
438
439
[2003]440   SUBROUTINE DDPDD( ydda, yddb )
441      !!----------------------------------------------------------------------
442      !!               ***  ROUTINE DDPDD ***
[3764]443      !!
[2003]444      !! ** Purpose : Add a scalar element to a sum
445      !!
[3764]446      !!
447      !! ** Method  : The code uses the compensated summation with doublet
[2003]448      !!              (sum,error) emulated useing complex numbers. ydda is the
[3764]449      !!               scalar to add to the summ yddb
[2003]450      !!
[3764]451      !! ** Action  : This does only work for MPI.
452      !!
[2003]453      !! References : Using Acurate Arithmetics to Improve Numerical
454      !!              Reproducibility and Sability in Parallel Applications
[3764]455      !!              Yun HE and Chris H. Q. DING, Journal of Supercomputing 18, 259-277, 2001
[2003]456      !!----------------------------------------------------------------------
[2307]457      COMPLEX(wp), INTENT(in   ) ::   ydda
458      COMPLEX(wp), INTENT(inout) ::   yddb
459      !
[2003]460      REAL(wp) :: zerr, zt1, zt2  ! local work variables
[2307]461      !!-----------------------------------------------------------------------
462      !
[2003]463      ! Compute ydda + yddb using Knuth's trick.
[2307]464      zt1  = REAL(ydda) + REAL(yddb)
465      zerr = zt1 - REAL(ydda)
466      zt2  = ( (REAL(yddb) - zerr) + (REAL(ydda) - (zt1 - zerr)) )   &
467         &   + AIMAG(ydda)         + AIMAG(yddb)
468      !
[2003]469      ! The result is t1 + t2, after normalization.
[2307]470      yddb = CMPLX( zt1 + zt2, zt2 - ((zt1 + zt2) - zt1), wp )
471      !
[2003]472   END SUBROUTINE DDPDD
473
474#if defined key_nosignedzero
[2307]475   !!----------------------------------------------------------------------
476   !!   'key_nosignedzero'                                         F90 SIGN
477   !!----------------------------------------------------------------------
[3764]478
[2307]479   FUNCTION SIGN_SCALAR( pa, pb )
[2003]480      !!-----------------------------------------------------------------------
481      !!                  ***  FUNCTION SIGN_SCALAR  ***
482      !!
483      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
484      !!-----------------------------------------------------------------------
485      REAL(wp) :: pa,pb          ! input
[2307]486      REAL(wp) :: SIGN_SCALAR    ! result
487      !!-----------------------------------------------------------------------
488      IF ( pb >= 0.e0) THEN   ;   SIGN_SCALAR = ABS(pa)
489      ELSE                    ;   SIGN_SCALAR =-ABS(pa)
[2003]490      ENDIF
491   END FUNCTION SIGN_SCALAR
492
[2307]493
[3764]494   FUNCTION SIGN_ARRAY_1D( pa, pb )
[2003]495      !!-----------------------------------------------------------------------
496      !!                  ***  FUNCTION SIGN_ARRAY_1D  ***
497      !!
498      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
499      !!-----------------------------------------------------------------------
[2307]500      REAL(wp) :: pa,pb(:)                   ! input
[2003]501      REAL(wp) :: SIGN_ARRAY_1D(SIZE(pb,1))  ! result
[2307]502      !!-----------------------------------------------------------------------
503      WHERE ( pb >= 0.e0 )   ;   SIGN_ARRAY_1D = ABS(pa)
504      ELSEWHERE              ;   SIGN_ARRAY_1D =-ABS(pa)
[2003]505      END WHERE
506   END FUNCTION SIGN_ARRAY_1D
507
[2307]508
[3764]509   FUNCTION SIGN_ARRAY_2D(pa,pb)
[2003]510      !!-----------------------------------------------------------------------
511      !!                  ***  FUNCTION SIGN_ARRAY_2D  ***
512      !!
513      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
514      !!-----------------------------------------------------------------------
515      REAL(wp) :: pa,pb(:,:)      ! input
516      REAL(wp) :: SIGN_ARRAY_2D(SIZE(pb,1),SIZE(pb,2))  ! result
[2307]517      !!-----------------------------------------------------------------------
518      WHERE ( pb >= 0.e0 )   ;   SIGN_ARRAY_2D = ABS(pa)
519      ELSEWHERE              ;   SIGN_ARRAY_2D =-ABS(pa)
[2003]520      END WHERE
521   END FUNCTION SIGN_ARRAY_2D
522
[3764]523   FUNCTION SIGN_ARRAY_3D(pa,pb)
[2003]524      !!-----------------------------------------------------------------------
525      !!                  ***  FUNCTION SIGN_ARRAY_3D  ***
526      !!
527      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
528      !!-----------------------------------------------------------------------
529      REAL(wp) :: pa,pb(:,:,:)      ! input
530      REAL(wp) :: SIGN_ARRAY_3D(SIZE(pb,1),SIZE(pb,2),SIZE(pb,3))  ! result
[2307]531      !!-----------------------------------------------------------------------
532      WHERE ( pb >= 0.e0 )   ;   SIGN_ARRAY_3D = ABS(pa)
533      ELSEWHERE              ;   SIGN_ARRAY_3D =-ABS(pa)
[2003]534      END WHERE
535   END FUNCTION SIGN_ARRAY_3D
536
[2307]537
[3764]538   FUNCTION SIGN_ARRAY_1D_A(pa,pb)
[2003]539      !!-----------------------------------------------------------------------
540      !!                  ***  FUNCTION SIGN_ARRAY_1D_A  ***
541      !!
542      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
543      !!-----------------------------------------------------------------------
544      REAL(wp) :: pa(:),pb(:)      ! input
[2307]545      REAL(wp) :: SIGN_ARRAY_1D_A(SIZE(pb,1))  ! result
546      !!-----------------------------------------------------------------------
547      WHERE ( pb >= 0.e0 )   ;   SIGN_ARRAY_1D_A = ABS(pa)
548      ELSEWHERE              ;   SIGN_ARRAY_1D_A =-ABS(pa)
[2003]549      END WHERE
550   END FUNCTION SIGN_ARRAY_1D_A
551
[2307]552
[3764]553   FUNCTION SIGN_ARRAY_2D_A(pa,pb)
[2003]554      !!-----------------------------------------------------------------------
555      !!                  ***  FUNCTION SIGN_ARRAY_2D_A  ***
556      !!
557      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
558      !!-----------------------------------------------------------------------
559      REAL(wp) :: pa(:,:),pb(:,:)      ! input
560      REAL(wp) :: SIGN_ARRAY_2D_A(SIZE(pb,1),SIZE(pb,2))  ! result
[2307]561      !!-----------------------------------------------------------------------
562      WHERE ( pb >= 0.e0 )   ;   SIGN_ARRAY_2D_A = ABS(pa)
563      ELSEWHERE              ;   SIGN_ARRAY_2D_A =-ABS(pa)
[2003]564      END WHERE
565   END FUNCTION SIGN_ARRAY_2D_A
566
[2307]567
[3764]568   FUNCTION SIGN_ARRAY_3D_A(pa,pb)
[2003]569      !!-----------------------------------------------------------------------
570      !!                  ***  FUNCTION SIGN_ARRAY_3D_A  ***
571      !!
572      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
573      !!-----------------------------------------------------------------------
574      REAL(wp) :: pa(:,:,:),pb(:,:,:)  ! input
575      REAL(wp) :: SIGN_ARRAY_3D_A(SIZE(pb,1),SIZE(pb,2),SIZE(pb,3)) ! result
[2307]576      !!-----------------------------------------------------------------------
577      WHERE ( pb >= 0.e0 )   ;   SIGN_ARRAY_3D_A = ABS(pa)
578      ELSEWHERE              ;   SIGN_ARRAY_3D_A =-ABS(pa)
[2003]579      END WHERE
580   END FUNCTION SIGN_ARRAY_3D_A
581
[2307]582
[3764]583   FUNCTION SIGN_ARRAY_1D_B(pa,pb)
[2003]584      !!-----------------------------------------------------------------------
585      !!                  ***  FUNCTION SIGN_ARRAY_1D_B  ***
586      !!
587      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
588      !!-----------------------------------------------------------------------
589      REAL(wp) :: pa(:),pb      ! input
590      REAL(wp) :: SIGN_ARRAY_1D_B(SIZE(pa,1))  ! result
[2307]591      !!-----------------------------------------------------------------------
592      IF( pb >= 0.e0 ) THEN   ;   SIGN_ARRAY_1D_B = ABS(pa)
593      ELSE                    ;   SIGN_ARRAY_1D_B =-ABS(pa)
[2003]594      ENDIF
595   END FUNCTION SIGN_ARRAY_1D_B
596
[2307]597
[3764]598   FUNCTION SIGN_ARRAY_2D_B(pa,pb)
[2003]599      !!-----------------------------------------------------------------------
600      !!                  ***  FUNCTION SIGN_ARRAY_2D_B  ***
601      !!
602      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
603      !!-----------------------------------------------------------------------
604      REAL(wp) :: pa(:,:),pb      ! input
605      REAL(wp) :: SIGN_ARRAY_2D_B(SIZE(pa,1),SIZE(pa,2))  ! result
[2307]606      !!-----------------------------------------------------------------------
607      IF( pb >= 0.e0 ) THEN   ;   SIGN_ARRAY_2D_B = ABS(pa)
608      ELSE                    ;   SIGN_ARRAY_2D_B =-ABS(pa)
[2003]609      ENDIF
610   END FUNCTION SIGN_ARRAY_2D_B
611
[2307]612
[3764]613   FUNCTION SIGN_ARRAY_3D_B(pa,pb)
[2003]614      !!-----------------------------------------------------------------------
615      !!                  ***  FUNCTION SIGN_ARRAY_3D_B  ***
616      !!
617      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
618      !!-----------------------------------------------------------------------
619      REAL(wp) :: pa(:,:,:),pb      ! input
620      REAL(wp) :: SIGN_ARRAY_3D_B(SIZE(pa,1),SIZE(pa,2),SIZE(pa,3))  ! result
[2307]621      !!-----------------------------------------------------------------------
622      IF( pb >= 0.e0 ) THEN   ;   SIGN_ARRAY_3D_B = ABS(pa)
623      ELSE                    ;   SIGN_ARRAY_3D_B =-ABS(pa)
[2003]624      ENDIF
625   END FUNCTION SIGN_ARRAY_3D_B
626#endif
627
[2307]628   !!======================================================================
[2003]629END MODULE lib_fortran
Note: See TracBrowser for help on using the repository browser.