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 on Ticket #936 – Attachment – NEMO

Ticket #936: lib_fortran.F90

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