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 branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC – NEMO

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/lib_fortran.F90 @ 4409

Last change on this file since 4409 was 3187, checked in by spickles2, 12 years ago

Stephen Pickles, 28 Nov 2011.
First commit of dCSE NEMO project work, part 1 - index re-ordering,
OPA_SRC top level only. Includes fix for sub-optimal auto-partitioning
in nemogcm.F90.

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