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

Last change on this file since 3421 was 3421, checked in by charris, 9 years ago

#936 Changes as suggested by Gurvan et al (testing details in the ticket). Note that for ORCA1 results will change with either nn_closea=0 or 1 because the Caspian Sea is now coded as a closed sea.

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