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

Last change on this file since 3294 was 3294, checked in by rblod, 12 years ago

Merge of 3.4beta into the trunk
  • Property svn:keywords set to Id
File size: 18.5 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   !!----------------------------------------------------------------------
8
9   !!----------------------------------------------------------------------
10   !!   glob_sum    : generic interface for global masked summation over
11   !!                 the interior domain for 1 or 2 2D or 3D arrays
12   !!                 it works only for T points   
13   !!   SIGN        : generic interface for SIGN to overwrite f95 behaviour
14   !!                 of intrinsinc sign function
15   !!----------------------------------------------------------------------
16   USE par_oce          ! Ocean parameter
17   USE lib_mpp          ! distributed memory computing
18   USE dom_oce          ! ocean domain
19   USE in_out_manager   ! I/O manager
20
21   IMPLICIT NONE
22   PRIVATE
23
24   PUBLIC glob_sum
25#if defined key_nosignedzero
26   PUBLIC SIGN
27#endif
28
29   INTERFACE glob_sum
30      MODULE PROCEDURE glob_sum_2d, glob_sum_3d,glob_sum_2d_a, glob_sum_3d_a 
31   END INTERFACE
32
33#if defined key_nosignedzero   
34   INTERFACE SIGN
35      MODULE PROCEDURE SIGN_SCALAR, SIGN_ARRAY_1D, SIGN_ARRAY_2D, SIGN_ARRAY_3D,   &
36         &             SIGN_ARRAY_1D_A, SIGN_ARRAY_2D_A, SIGN_ARRAY_3D_A,          & 
37         &             SIGN_ARRAY_1D_B, SIGN_ARRAY_2D_B, SIGN_ARRAY_3D_B 
38   END INTERFACE
39#endif
40
41   !!----------------------------------------------------------------------
42   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
43   !! $Id$
44   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
45   !!----------------------------------------------------------------------
46CONTAINS 
47
48#if ! defined key_mpp_rep
49   FUNCTION glob_sum_2d( ptab ) 
50      !!-----------------------------------------------------------------------
51      !!                  ***  FUNCTION  glob_sum_2D  ***
52      !!
53      !! ** Purpose : perform a masked sum on the inner global domain of a 2D array
54      !!-----------------------------------------------------------------------
55      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptab          ! input 2D array
56      REAL(wp)                             ::   glob_sum_2d   ! global masked sum
57      !!-----------------------------------------------------------------------
58      !
59      glob_sum_2d = SUM( ptab(:,:)*tmask_i(:,:) )
60      IF( lk_mpp )   CALL mpp_sum( glob_sum_2d )
61      !
62   END FUNCTION glob_sum_2d
63   
64   
65   FUNCTION glob_sum_3d( ptab ) 
66      !!-----------------------------------------------------------------------
67      !!                  ***  FUNCTION  glob_sum_3D  ***
68      !!
69      !! ** Purpose : perform a masked sum on the inner global domain of a 3D array
70      !!-----------------------------------------------------------------------
71      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   ptab          ! input 3D array
72      REAL(wp)                               ::   glob_sum_3d   ! global masked sum
73      !!
74      INTEGER :: jk
75      !!-----------------------------------------------------------------------
76      !
77      glob_sum_3d = 0.e0
78      DO jk = 1, jpk
79         glob_sum_3d = glob_sum_3d + SUM( ptab(:,:,jk)*tmask_i(:,:) )
80      END DO
81      IF( lk_mpp )   CALL mpp_sum( glob_sum_3d )
82      !
83   END FUNCTION glob_sum_3d
84
85
86   FUNCTION glob_sum_2d_a( ptab1, ptab2 ) 
87      !!-----------------------------------------------------------------------
88      !!                  ***  FUNCTION  glob_sum_2D _a ***
89      !!
90      !! ** Purpose : perform a masked sum on the inner global domain of two 2D array
91      !!-----------------------------------------------------------------------
92      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptab1, ptab2    ! input 2D array
93      REAL(wp)            , DIMENSION(2)   ::   glob_sum_2d_a   ! global masked sum
94      !!-----------------------------------------------------------------------
95      !             
96      glob_sum_2d_a(1) = SUM( ptab1(:,:)*tmask_i(:,:) )
97      glob_sum_2d_a(2) = SUM( ptab2(:,:)*tmask_i(:,:) )
98      IF( lk_mpp )   CALL mpp_sum( glob_sum_2d_a, 2 )
99      !
100   END FUNCTION glob_sum_2d_a
101 
102 
103   FUNCTION glob_sum_3d_a( ptab1, ptab2 ) 
104      !!-----------------------------------------------------------------------
105      !!                  ***  FUNCTION  glob_sum_3D_a ***
106      !!
107      !! ** Purpose : perform a masked sum on the inner global domain of two 3D array
108      !!-----------------------------------------------------------------------
109      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   ptab1, ptab2    ! input 3D array
110      REAL(wp)            , DIMENSION(2)     ::   glob_sum_3d_a   ! global masked sum
111      !!
112      INTEGER :: jk
113      !!-----------------------------------------------------------------------
114      !
115      glob_sum_3d_a(:) = 0.e0
116      DO jk = 1, jpk
117         glob_sum_3d_a(1) = glob_sum_3d_a(1) + SUM( ptab1(:,:,jk)*tmask_i(:,:) )
118         glob_sum_3d_a(2) = glob_sum_3d_a(2) + SUM( ptab2(:,:,jk)*tmask_i(:,:) )
119      END DO
120      IF( lk_mpp )   CALL mpp_sum( glob_sum_3d_a, 2 )
121      !
122   END FUNCTION glob_sum_3d_a
123
124#else 
125   !!----------------------------------------------------------------------
126   !!   'key_mpp_rep'                                   MPP reproducibility
127   !!----------------------------------------------------------------------
128   
129   FUNCTION glob_sum_2d( ptab ) 
130      !!----------------------------------------------------------------------
131      !!                  ***  FUNCTION  glob_sum_2d ***
132      !!
133      !! ** Purpose : perform a sum in calling DDPDD routine
134      !!----------------------------------------------------------------------
135      REAL(wp), INTENT(in), DIMENSION(jpi,jpj) ::   ptab
136      REAL(wp)                                 ::   glob_sum_2d   ! global masked sum
137      !!
138      COMPLEX(wp)::   ctmp
139      REAL(wp)   ::   ztmp
140      INTEGER    ::   ji, jj   ! dummy loop indices
141      !!-----------------------------------------------------------------------
142      !
143      ztmp = 0.e0
144      ctmp = CMPLX( 0.e0, 0.e0, wp )
145      DO jj = 1, jpj
146         DO ji =1, jpi
147         ztmp =  ptab(ji,jj) * tmask_i(ji,jj)
148         CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
149         END DO
150      END DO
151      IF( lk_mpp )   CALL mpp_sum( ctmp )   ! sum over the global domain
152      glob_sum_2d = REAL(ctmp,wp)
153      !
154   END FUNCTION glob_sum_2d   
155
156
157   FUNCTION glob_sum_3d( ptab ) 
158      !!----------------------------------------------------------------------
159      !!                  ***  FUNCTION  glob_sum_3d ***
160      !!
161      !! ** Purpose : perform a sum on a 3D array in calling DDPDD routine
162      !!----------------------------------------------------------------------
163      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   ptab
164      REAL(wp)                                     ::   glob_sum_3d   ! global masked sum
165      !!
166      COMPLEX(wp)::   ctmp
167      REAL(wp)   ::   ztmp
168      INTEGER    ::   ji, jj, jk   ! dummy loop indices
169      !!-----------------------------------------------------------------------
170      !
171      ztmp = 0.e0
172      ctmp = CMPLX( 0.e0, 0.e0, wp )
173      DO jk = 1, jpk
174         DO jj = 1, jpj
175            DO ji =1, jpi
176            ztmp =  ptab(ji,jj,jk) * tmask_i(ji,jj)
177            CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
178            END DO
179         END DO   
180      END DO
181      IF( lk_mpp )   CALL mpp_sum( ctmp )   ! sum over the global domain
182      glob_sum_3d = REAL(ctmp,wp)
183      !
184   END FUNCTION glob_sum_3d   
185
186
187   FUNCTION glob_sum_2d_a( ptab1, ptab2 ) 
188      !!----------------------------------------------------------------------
189      !!                  ***  FUNCTION  glob_sum_2d_a ***
190      !!
191      !! ** Purpose : perform a sum on two 2D arrays in calling DDPDD routine
192      !!----------------------------------------------------------------------
193      REAL(wp), INTENT(in), DIMENSION(jpi,jpj) ::   ptab1, ptab2
194      REAL(wp)                                 ::   glob_sum_2d_a   ! global masked sum
195      !!
196      COMPLEX(wp)::   ctmp
197      REAL(wp)   ::   ztmp
198      INTEGER    ::   ji, jj   ! dummy loop indices
199      !!-----------------------------------------------------------------------
200      !
201      ztmp = 0.e0
202      ctmp = CMPLX( 0.e0, 0.e0, wp )
203      DO jj = 1, jpj
204         DO ji =1, jpi
205         ztmp =  ptab1(ji,jj) * tmask_i(ji,jj)
206         CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
207         ztmp =  ptab2(ji,jj) * tmask_i(ji,jj)
208         CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
209         END DO
210      END DO
211      IF( lk_mpp )   CALL mpp_sum( ctmp )   ! sum over the global domain
212      glob_sum_2d_a = REAL(ctmp,wp)
213      !
214   END FUNCTION glob_sum_2d_a   
215
216
217   FUNCTION glob_sum_3d_a( ptab1, ptab2 ) 
218      !!----------------------------------------------------------------------
219      !!                  ***  FUNCTION  glob_sum_3d_a ***
220      !!
221      !! ** Purpose : perform a sum on two 3D array in calling DDPDD routine
222      !!----------------------------------------------------------------------
223      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   ptab1, ptab2
224      REAL(wp)                                     ::   glob_sum_3d_a   ! global masked sum
225      !!
226      COMPLEX(wp)::   ctmp
227      REAL(wp)   ::   ztmp
228      INTEGER    ::   ji, jj, jk   ! dummy loop indices
229      !!-----------------------------------------------------------------------
230      !
231      ztmp = 0.e0
232      ctmp = CMPLX( 0.e0, 0.e0, wp )
233      DO jk = 1, jpk
234         DO jj = 1, jpj
235            DO ji =1, jpi
236            ztmp =  ptab1(ji,jj,jk) * tmask_i(ji,jj)
237            CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
238            ztmp =  ptab2(ji,jj,jk) * tmask_i(ji,jj)
239            CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
240            END DO
241         END DO   
242      END DO
243      IF( lk_mpp )   CALL mpp_sum( ctmp )   ! sum over the global domain
244      glob_sum_3d_a = REAL(ctmp,wp)
245      !
246   END FUNCTION glob_sum_3d_a   
247
248
249   SUBROUTINE DDPDD( ydda, yddb )
250      !!----------------------------------------------------------------------
251      !!               ***  ROUTINE DDPDD ***
252      !!         
253      !! ** Purpose : Add a scalar element to a sum
254      !!             
255      !!
256      !! ** Method  : The code uses the compensated summation with doublet
257      !!              (sum,error) emulated useing complex numbers. ydda is the
258      !!               scalar to add to the summ yddb
259      !!
260      !! ** Action  : This does only work for MPI.
261      !!
262      !! References : Using Acurate Arithmetics to Improve Numerical
263      !!              Reproducibility and Sability in Parallel Applications
264      !!              Yun HE and Chris H. Q. DING, Journal of Supercomputing 18, 259-277, 2001
265      !!----------------------------------------------------------------------
266      COMPLEX(wp), INTENT(in   ) ::   ydda
267      COMPLEX(wp), INTENT(inout) ::   yddb
268      !
269      REAL(wp) :: zerr, zt1, zt2  ! local work variables
270      !!-----------------------------------------------------------------------
271      !
272      ! Compute ydda + yddb using Knuth's trick.
273      zt1  = REAL(ydda) + REAL(yddb)
274      zerr = zt1 - REAL(ydda)
275      zt2  = ( (REAL(yddb) - zerr) + (REAL(ydda) - (zt1 - zerr)) )   &
276         &   + AIMAG(ydda)         + AIMAG(yddb)
277      !
278      ! The result is t1 + t2, after normalization.
279      yddb = CMPLX( zt1 + zt2, zt2 - ((zt1 + zt2) - zt1), wp )
280      !
281   END SUBROUTINE DDPDD
282#endif
283
284#if defined key_nosignedzero
285   !!----------------------------------------------------------------------
286   !!   'key_nosignedzero'                                         F90 SIGN
287   !!----------------------------------------------------------------------
288   
289   FUNCTION SIGN_SCALAR( pa, pb )
290      !!-----------------------------------------------------------------------
291      !!                  ***  FUNCTION SIGN_SCALAR  ***
292      !!
293      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
294      !!-----------------------------------------------------------------------
295      REAL(wp) :: pa,pb          ! input
296      REAL(wp) :: SIGN_SCALAR    ! result
297      !!-----------------------------------------------------------------------
298      IF ( pb >= 0.e0) THEN   ;   SIGN_SCALAR = ABS(pa)
299      ELSE                    ;   SIGN_SCALAR =-ABS(pa)
300      ENDIF
301   END FUNCTION SIGN_SCALAR
302
303
304   FUNCTION SIGN_ARRAY_1D( pa, pb ) 
305      !!-----------------------------------------------------------------------
306      !!                  ***  FUNCTION SIGN_ARRAY_1D  ***
307      !!
308      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
309      !!-----------------------------------------------------------------------
310      REAL(wp) :: pa,pb(:)                   ! input
311      REAL(wp) :: SIGN_ARRAY_1D(SIZE(pb,1))  ! result
312      !!-----------------------------------------------------------------------
313      WHERE ( pb >= 0.e0 )   ;   SIGN_ARRAY_1D = ABS(pa)
314      ELSEWHERE              ;   SIGN_ARRAY_1D =-ABS(pa)
315      END WHERE
316   END FUNCTION SIGN_ARRAY_1D
317
318
319   FUNCTION SIGN_ARRAY_2D(pa,pb) 
320      !!-----------------------------------------------------------------------
321      !!                  ***  FUNCTION SIGN_ARRAY_2D  ***
322      !!
323      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
324      !!-----------------------------------------------------------------------
325      REAL(wp) :: pa,pb(:,:)      ! input
326      REAL(wp) :: SIGN_ARRAY_2D(SIZE(pb,1),SIZE(pb,2))  ! result
327      !!-----------------------------------------------------------------------
328      WHERE ( pb >= 0.e0 )   ;   SIGN_ARRAY_2D = ABS(pa)
329      ELSEWHERE              ;   SIGN_ARRAY_2D =-ABS(pa)
330      END WHERE
331   END FUNCTION SIGN_ARRAY_2D
332
333   FUNCTION SIGN_ARRAY_3D(pa,pb) 
334      !!-----------------------------------------------------------------------
335      !!                  ***  FUNCTION SIGN_ARRAY_3D  ***
336      !!
337      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
338      !!-----------------------------------------------------------------------
339      REAL(wp) :: pa,pb(:,:,:)      ! input
340      REAL(wp) :: SIGN_ARRAY_3D(SIZE(pb,1),SIZE(pb,2),SIZE(pb,3))  ! result
341      !!-----------------------------------------------------------------------
342      WHERE ( pb >= 0.e0 )   ;   SIGN_ARRAY_3D = ABS(pa)
343      ELSEWHERE              ;   SIGN_ARRAY_3D =-ABS(pa)
344      END WHERE
345   END FUNCTION SIGN_ARRAY_3D
346
347
348   FUNCTION SIGN_ARRAY_1D_A(pa,pb) 
349      !!-----------------------------------------------------------------------
350      !!                  ***  FUNCTION SIGN_ARRAY_1D_A  ***
351      !!
352      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
353      !!-----------------------------------------------------------------------
354      REAL(wp) :: pa(:),pb(:)      ! input
355      REAL(wp) :: SIGN_ARRAY_1D_A(SIZE(pb,1))  ! result
356      !!-----------------------------------------------------------------------
357      WHERE ( pb >= 0.e0 )   ;   SIGN_ARRAY_1D_A = ABS(pa)
358      ELSEWHERE              ;   SIGN_ARRAY_1D_A =-ABS(pa)
359      END WHERE
360   END FUNCTION SIGN_ARRAY_1D_A
361
362
363   FUNCTION SIGN_ARRAY_2D_A(pa,pb) 
364      !!-----------------------------------------------------------------------
365      !!                  ***  FUNCTION SIGN_ARRAY_2D_A  ***
366      !!
367      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
368      !!-----------------------------------------------------------------------
369      REAL(wp) :: pa(:,:),pb(:,:)      ! input
370      REAL(wp) :: SIGN_ARRAY_2D_A(SIZE(pb,1),SIZE(pb,2))  ! result
371      !!-----------------------------------------------------------------------
372      WHERE ( pb >= 0.e0 )   ;   SIGN_ARRAY_2D_A = ABS(pa)
373      ELSEWHERE              ;   SIGN_ARRAY_2D_A =-ABS(pa)
374      END WHERE
375   END FUNCTION SIGN_ARRAY_2D_A
376
377
378   FUNCTION SIGN_ARRAY_3D_A(pa,pb) 
379      !!-----------------------------------------------------------------------
380      !!                  ***  FUNCTION SIGN_ARRAY_3D_A  ***
381      !!
382      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
383      !!-----------------------------------------------------------------------
384      REAL(wp) :: pa(:,:,:),pb(:,:,:)  ! input
385      REAL(wp) :: SIGN_ARRAY_3D_A(SIZE(pb,1),SIZE(pb,2),SIZE(pb,3)) ! result
386      !!-----------------------------------------------------------------------
387      WHERE ( pb >= 0.e0 )   ;   SIGN_ARRAY_3D_A = ABS(pa)
388      ELSEWHERE              ;   SIGN_ARRAY_3D_A =-ABS(pa)
389      END WHERE
390   END FUNCTION SIGN_ARRAY_3D_A
391
392
393   FUNCTION SIGN_ARRAY_1D_B(pa,pb) 
394      !!-----------------------------------------------------------------------
395      !!                  ***  FUNCTION SIGN_ARRAY_1D_B  ***
396      !!
397      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
398      !!-----------------------------------------------------------------------
399      REAL(wp) :: pa(:),pb      ! input
400      REAL(wp) :: SIGN_ARRAY_1D_B(SIZE(pa,1))  ! result
401      !!-----------------------------------------------------------------------
402      IF( pb >= 0.e0 ) THEN   ;   SIGN_ARRAY_1D_B = ABS(pa)
403      ELSE                    ;   SIGN_ARRAY_1D_B =-ABS(pa)
404      ENDIF
405   END FUNCTION SIGN_ARRAY_1D_B
406
407
408   FUNCTION SIGN_ARRAY_2D_B(pa,pb) 
409      !!-----------------------------------------------------------------------
410      !!                  ***  FUNCTION SIGN_ARRAY_2D_B  ***
411      !!
412      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
413      !!-----------------------------------------------------------------------
414      REAL(wp) :: pa(:,:),pb      ! input
415      REAL(wp) :: SIGN_ARRAY_2D_B(SIZE(pa,1),SIZE(pa,2))  ! result
416      !!-----------------------------------------------------------------------
417      IF( pb >= 0.e0 ) THEN   ;   SIGN_ARRAY_2D_B = ABS(pa)
418      ELSE                    ;   SIGN_ARRAY_2D_B =-ABS(pa)
419      ENDIF
420   END FUNCTION SIGN_ARRAY_2D_B
421
422
423   FUNCTION SIGN_ARRAY_3D_B(pa,pb) 
424      !!-----------------------------------------------------------------------
425      !!                  ***  FUNCTION SIGN_ARRAY_3D_B  ***
426      !!
427      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
428      !!-----------------------------------------------------------------------
429      REAL(wp) :: pa(:,:,:),pb      ! input
430      REAL(wp) :: SIGN_ARRAY_3D_B(SIZE(pa,1),SIZE(pa,2),SIZE(pa,3))  ! result
431      !!-----------------------------------------------------------------------
432      IF( pb >= 0.e0 ) THEN   ;   SIGN_ARRAY_3D_B = ABS(pa)
433      ELSE                    ;   SIGN_ARRAY_3D_B =-ABS(pa)
434      ENDIF
435   END FUNCTION SIGN_ARRAY_3D_B
436#endif
437
438   !!======================================================================
439END MODULE lib_fortran
Note: See TracBrowser for help on using the repository browser.