source: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/lib_fortran.F90 @ 10397

Last change on this file since 10397 was 10397, checked in by smasson, 23 months ago

dev_r10164_HPC09_ESIWACE_PREP_MERGE: action 5 introduce mpp_delay_sum, see #2133

  • Property svn:keywords set to Id
File size: 15.2 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   !!            3.4  !  2013-06  (C. Rousset)  add glob_min, glob_max
8   !!                                           + 3d dim. of input is fexible (jpk, jpl...)
9   !!            4.0  !  2016-06  (T. Lovato)  double precision global sum by default
10   !!----------------------------------------------------------------------
11
12   !!----------------------------------------------------------------------
13   !!   glob_sum    : generic interface for global masked summation over
14   !!                 the interior domain for 1 or 2 2D or 3D arrays
15   !!                 it works only for T points
16   !!   SIGN        : generic interface for SIGN to overwrite f95 behaviour
17   !!                 of intrinsinc sign function
18   !!----------------------------------------------------------------------
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
23
24   IMPLICIT NONE
25   PRIVATE
26
27   PUBLIC   glob_sum      ! used in many places (masked with tmask_i)
28   PUBLIC   glob_sum_full ! used in many places (masked with tmask_h, ie only over the halos)
29   PUBLIC   local_sum     ! used in trcrad, local operation before glob_sum_delay
30   PUBLIC   DDPDD         ! also used in closea module
31   PUBLIC   glob_min, glob_max
32#if defined key_nosignedzero
33   PUBLIC SIGN
34#endif
35
36   INTERFACE glob_sum
37      MODULE PROCEDURE glob_sum_1d, glob_sum_2d, glob_sum_3d
38   END INTERFACE
39   INTERFACE glob_sum_full
40      MODULE PROCEDURE glob_sum_full_2d, glob_sum_full_3d
41   END INTERFACE
42   INTERFACE local_sum
43      MODULE PROCEDURE local_sum_2d, local_sum_3d
44   END INTERFACE
45   INTERFACE glob_min
46      MODULE PROCEDURE glob_min_2d, glob_min_3d
47   END INTERFACE
48   INTERFACE glob_max
49      MODULE PROCEDURE glob_max_2d, glob_max_3d
50   END INTERFACE
51
52#if defined key_nosignedzero
53   INTERFACE SIGN
54      MODULE PROCEDURE SIGN_SCALAR, SIGN_ARRAY_1D, SIGN_ARRAY_2D, SIGN_ARRAY_3D,   &
55         &             SIGN_ARRAY_1D_A, SIGN_ARRAY_2D_A, SIGN_ARRAY_3D_A,          &
56         &             SIGN_ARRAY_1D_B, SIGN_ARRAY_2D_B, SIGN_ARRAY_3D_B
57   END INTERFACE
58#endif
59
60   !!----------------------------------------------------------------------
61   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
62   !! $Id$
63   !! Software governed by the CeCILL license (see ./LICENSE)
64   !!----------------------------------------------------------------------
65CONTAINS
66
67#  define GLOBSUM_CODE
68
69#     define DIM_1d
70#     define FUNCTION_GLOBSUM           glob_sum_1d
71#     include "lib_fortran_generic.h90"
72#     undef FUNCTION_GLOBSUM
73#     undef DIM_1d
74
75#     define DIM_2d
76#     define OPERATION_GLOBSUM
77#     define FUNCTION_GLOBSUM           glob_sum_2d
78#     include "lib_fortran_generic.h90"
79#     undef FUNCTION_GLOBSUM
80#     undef OPERATION_GLOBSUM
81#     define OPERATION_FULL_GLOBSUM
82#     define FUNCTION_GLOBSUM           glob_sum_full_2d
83#     include "lib_fortran_generic.h90"
84#     undef FUNCTION_GLOBSUM
85#     undef OPERATION_FULL_GLOBSUM
86#     undef DIM_2d
87
88#     define DIM_3d
89#     define OPERATION_GLOBSUM
90#     define FUNCTION_GLOBSUM           glob_sum_3d
91#     include "lib_fortran_generic.h90"
92#     undef FUNCTION_GLOBSUM
93#     undef OPERATION_GLOBSUM
94#     define OPERATION_FULL_GLOBSUM
95#     define FUNCTION_GLOBSUM           glob_sum_full_3d
96#     include "lib_fortran_generic.h90"
97#     undef FUNCTION_GLOBSUM
98#     undef OPERATION_FULL_GLOBSUM
99#     undef DIM_3d
100
101#  undef GLOBSUM_CODE
102
103
104#  define GLOBMINMAX_CODE
105
106#     define DIM_2d
107#     define OPERATION_GLOBMIN
108#     define FUNCTION_GLOBMINMAX           glob_min_2d
109#     include "lib_fortran_generic.h90"
110#     undef FUNCTION_GLOBMINMAX
111#     undef OPERATION_GLOBMIN
112#     define OPERATION_GLOBMAX
113#     define FUNCTION_GLOBMINMAX           glob_max_2d
114#     include "lib_fortran_generic.h90"
115#     undef FUNCTION_GLOBMINMAX
116#     undef OPERATION_GLOBMAX
117#     undef DIM_2d
118
119#     define DIM_3d
120#     define OPERATION_GLOBMIN
121#     define FUNCTION_GLOBMINMAX           glob_min_3d
122#     include "lib_fortran_generic.h90"
123#     undef FUNCTION_GLOBMINMAX
124#     undef OPERATION_GLOBMIN
125#     define OPERATION_GLOBMAX
126#     define FUNCTION_GLOBMINMAX           glob_max_3d
127#     include "lib_fortran_generic.h90"
128#     undef FUNCTION_GLOBMINMAX
129#     undef OPERATION_GLOBMAX
130#     undef DIM_3d
131#  undef GLOBMINMAX_CODE
132
133!                          ! FUNCTION local_sum !
134
135   FUNCTION local_sum_2d( ptab )
136      !!----------------------------------------------------------------------
137      REAL(wp),  INTENT(in   ) ::   ptab(:,:) ! array on which operation is applied
138      COMPLEX(wp)              ::  local_sum_2d
139      !
140      !!-----------------------------------------------------------------------
141      !
142      COMPLEX(wp)::   ctmp
143      REAL(wp)   ::   ztmp
144      INTEGER    ::   ji, jj    ! dummy loop indices
145      INTEGER    ::   ipi, ipj  ! dimensions
146      !!-----------------------------------------------------------------------
147      !
148      ipi = SIZE(ptab,1)   ! 1st dimension
149      ipj = SIZE(ptab,2)   ! 2nd dimension
150      !
151      ctmp = CMPLX( 0.e0, 0.e0, wp )   ! warning ctmp is cumulated
152
153      DO jj = 1, ipj
154         DO ji = 1, ipi
155            ztmp =  ptab(ji,jj) * tmask_i(ji,jj)
156            CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
157         END DO
158      END DO
159      !
160      local_sum_2d = ctmp
161       
162   END FUNCTION local_sum_2d
163
164   FUNCTION local_sum_3d( ptab )
165      !!----------------------------------------------------------------------
166      REAL(wp),  INTENT(in   ) ::   ptab(:,:,:) ! array on which operation is applied
167      COMPLEX(wp)              ::  local_sum_3d
168      !
169      !!-----------------------------------------------------------------------
170      !
171      COMPLEX(wp)::   ctmp
172      REAL(wp)   ::   ztmp
173      INTEGER    ::   ji, jj, jk   ! dummy loop indices
174      INTEGER    ::   ipi, ipj, ipk    ! dimensions
175      !!-----------------------------------------------------------------------
176      !
177      ipi = SIZE(ptab,1)   ! 1st dimension
178      ipj = SIZE(ptab,2)   ! 2nd dimension
179      ipk = SIZE(ptab,3)   ! 3rd dimension
180      !
181      ctmp = CMPLX( 0.e0, 0.e0, wp )   ! warning ctmp is cumulated
182
183      DO jk = 1, ipk
184        DO jj = 1, ipj
185          DO ji = 1, ipi
186             ztmp =  ptab(ji,jj,jk) * tmask_i(ji,jj)
187             CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
188          END DO
189        END DO
190      END DO
191      !
192      local_sum_3d = ctmp
193       
194   END FUNCTION local_sum_3d
195
196
197   SUBROUTINE DDPDD( ydda, yddb )
198      !!----------------------------------------------------------------------
199      !!               ***  ROUTINE DDPDD ***
200      !!
201      !! ** Purpose : Add a scalar element to a sum
202      !!
203      !!
204      !! ** Method  : The code uses the compensated summation with doublet
205      !!              (sum,error) emulated useing complex numbers. ydda is the
206      !!               scalar to add to the summ yddb
207      !!
208      !! ** Action  : This does only work for MPI.
209      !!
210      !! References : Using Acurate Arithmetics to Improve Numerical
211      !!              Reproducibility and Sability in Parallel Applications
212      !!              Yun HE and Chris H. Q. DING, Journal of Supercomputing 18, 259-277, 2001
213      !!----------------------------------------------------------------------
214      COMPLEX(wp), INTENT(in   ) ::   ydda
215      COMPLEX(wp), INTENT(inout) ::   yddb
216      !
217      REAL(wp) :: zerr, zt1, zt2  ! local work variables
218      !!-----------------------------------------------------------------------
219      !
220      ! Compute ydda + yddb using Knuth's trick.
221      zt1  = REAL(ydda) + REAL(yddb)
222      zerr = zt1 - REAL(ydda)
223      zt2  = ( (REAL(yddb) - zerr) + (REAL(ydda) - (zt1 - zerr)) )   &
224         &   + AIMAG(ydda)         + AIMAG(yddb)
225      !
226      ! The result is t1 + t2, after normalization.
227      yddb = CMPLX( zt1 + zt2, zt2 - ((zt1 + zt2) - zt1), wp )
228      !
229   END SUBROUTINE DDPDD
230
231#if defined key_nosignedzero
232   !!----------------------------------------------------------------------
233   !!   'key_nosignedzero'                                         F90 SIGN
234   !!----------------------------------------------------------------------
235
236   FUNCTION SIGN_SCALAR( pa, pb )
237      !!-----------------------------------------------------------------------
238      !!                  ***  FUNCTION SIGN_SCALAR  ***
239      !!
240      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
241      !!-----------------------------------------------------------------------
242      REAL(wp) :: pa,pb          ! input
243      REAL(wp) :: SIGN_SCALAR    ! result
244      !!-----------------------------------------------------------------------
245      IF ( pb >= 0.e0) THEN   ;   SIGN_SCALAR = ABS(pa)
246      ELSE                    ;   SIGN_SCALAR =-ABS(pa)
247      ENDIF
248   END FUNCTION SIGN_SCALAR
249
250
251   FUNCTION SIGN_ARRAY_1D( pa, pb )
252      !!-----------------------------------------------------------------------
253      !!                  ***  FUNCTION SIGN_ARRAY_1D  ***
254      !!
255      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
256      !!-----------------------------------------------------------------------
257      REAL(wp) :: pa,pb(:)                   ! input
258      REAL(wp) :: SIGN_ARRAY_1D(SIZE(pb,1))  ! result
259      !!-----------------------------------------------------------------------
260      WHERE ( pb >= 0.e0 )   ;   SIGN_ARRAY_1D = ABS(pa)
261      ELSEWHERE              ;   SIGN_ARRAY_1D =-ABS(pa)
262      END WHERE
263   END FUNCTION SIGN_ARRAY_1D
264
265
266   FUNCTION SIGN_ARRAY_2D(pa,pb)
267      !!-----------------------------------------------------------------------
268      !!                  ***  FUNCTION SIGN_ARRAY_2D  ***
269      !!
270      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
271      !!-----------------------------------------------------------------------
272      REAL(wp) :: pa,pb(:,:)      ! input
273      REAL(wp) :: SIGN_ARRAY_2D(SIZE(pb,1),SIZE(pb,2))  ! result
274      !!-----------------------------------------------------------------------
275      WHERE ( pb >= 0.e0 )   ;   SIGN_ARRAY_2D = ABS(pa)
276      ELSEWHERE              ;   SIGN_ARRAY_2D =-ABS(pa)
277      END WHERE
278   END FUNCTION SIGN_ARRAY_2D
279
280   FUNCTION SIGN_ARRAY_3D(pa,pb)
281      !!-----------------------------------------------------------------------
282      !!                  ***  FUNCTION SIGN_ARRAY_3D  ***
283      !!
284      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
285      !!-----------------------------------------------------------------------
286      REAL(wp) :: pa,pb(:,:,:)      ! input
287      REAL(wp) :: SIGN_ARRAY_3D(SIZE(pb,1),SIZE(pb,2),SIZE(pb,3))  ! result
288      !!-----------------------------------------------------------------------
289      WHERE ( pb >= 0.e0 )   ;   SIGN_ARRAY_3D = ABS(pa)
290      ELSEWHERE              ;   SIGN_ARRAY_3D =-ABS(pa)
291      END WHERE
292   END FUNCTION SIGN_ARRAY_3D
293
294
295   FUNCTION SIGN_ARRAY_1D_A(pa,pb)
296      !!-----------------------------------------------------------------------
297      !!                  ***  FUNCTION SIGN_ARRAY_1D_A  ***
298      !!
299      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
300      !!-----------------------------------------------------------------------
301      REAL(wp) :: pa(:),pb(:)      ! input
302      REAL(wp) :: SIGN_ARRAY_1D_A(SIZE(pb,1))  ! result
303      !!-----------------------------------------------------------------------
304      WHERE ( pb >= 0.e0 )   ;   SIGN_ARRAY_1D_A = ABS(pa)
305      ELSEWHERE              ;   SIGN_ARRAY_1D_A =-ABS(pa)
306      END WHERE
307   END FUNCTION SIGN_ARRAY_1D_A
308
309
310   FUNCTION SIGN_ARRAY_2D_A(pa,pb)
311      !!-----------------------------------------------------------------------
312      !!                  ***  FUNCTION SIGN_ARRAY_2D_A  ***
313      !!
314      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
315      !!-----------------------------------------------------------------------
316      REAL(wp) :: pa(:,:),pb(:,:)      ! input
317      REAL(wp) :: SIGN_ARRAY_2D_A(SIZE(pb,1),SIZE(pb,2))  ! result
318      !!-----------------------------------------------------------------------
319      WHERE ( pb >= 0.e0 )   ;   SIGN_ARRAY_2D_A = ABS(pa)
320      ELSEWHERE              ;   SIGN_ARRAY_2D_A =-ABS(pa)
321      END WHERE
322   END FUNCTION SIGN_ARRAY_2D_A
323
324
325   FUNCTION SIGN_ARRAY_3D_A(pa,pb)
326      !!-----------------------------------------------------------------------
327      !!                  ***  FUNCTION SIGN_ARRAY_3D_A  ***
328      !!
329      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
330      !!-----------------------------------------------------------------------
331      REAL(wp) :: pa(:,:,:),pb(:,:,:)  ! input
332      REAL(wp) :: SIGN_ARRAY_3D_A(SIZE(pb,1),SIZE(pb,2),SIZE(pb,3)) ! result
333      !!-----------------------------------------------------------------------
334      WHERE ( pb >= 0.e0 )   ;   SIGN_ARRAY_3D_A = ABS(pa)
335      ELSEWHERE              ;   SIGN_ARRAY_3D_A =-ABS(pa)
336      END WHERE
337   END FUNCTION SIGN_ARRAY_3D_A
338
339
340   FUNCTION SIGN_ARRAY_1D_B(pa,pb)
341      !!-----------------------------------------------------------------------
342      !!                  ***  FUNCTION SIGN_ARRAY_1D_B  ***
343      !!
344      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
345      !!-----------------------------------------------------------------------
346      REAL(wp) :: pa(:),pb      ! input
347      REAL(wp) :: SIGN_ARRAY_1D_B(SIZE(pa,1))  ! result
348      !!-----------------------------------------------------------------------
349      IF( pb >= 0.e0 ) THEN   ;   SIGN_ARRAY_1D_B = ABS(pa)
350      ELSE                    ;   SIGN_ARRAY_1D_B =-ABS(pa)
351      ENDIF
352   END FUNCTION SIGN_ARRAY_1D_B
353
354
355   FUNCTION SIGN_ARRAY_2D_B(pa,pb)
356      !!-----------------------------------------------------------------------
357      !!                  ***  FUNCTION SIGN_ARRAY_2D_B  ***
358      !!
359      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
360      !!-----------------------------------------------------------------------
361      REAL(wp) :: pa(:,:),pb      ! input
362      REAL(wp) :: SIGN_ARRAY_2D_B(SIZE(pa,1),SIZE(pa,2))  ! result
363      !!-----------------------------------------------------------------------
364      IF( pb >= 0.e0 ) THEN   ;   SIGN_ARRAY_2D_B = ABS(pa)
365      ELSE                    ;   SIGN_ARRAY_2D_B =-ABS(pa)
366      ENDIF
367   END FUNCTION SIGN_ARRAY_2D_B
368
369
370   FUNCTION SIGN_ARRAY_3D_B(pa,pb)
371      !!-----------------------------------------------------------------------
372      !!                  ***  FUNCTION SIGN_ARRAY_3D_B  ***
373      !!
374      !! ** Purpose : overwrite f95 behaviour of intrinsinc sign function
375      !!-----------------------------------------------------------------------
376      REAL(wp) :: pa(:,:,:),pb      ! input
377      REAL(wp) :: SIGN_ARRAY_3D_B(SIZE(pa,1),SIZE(pa,2),SIZE(pa,3))  ! result
378      !!-----------------------------------------------------------------------
379      IF( pb >= 0.e0 ) THEN   ;   SIGN_ARRAY_3D_B = ABS(pa)
380      ELSE                    ;   SIGN_ARRAY_3D_B =-ABS(pa)
381      ENDIF
382   END FUNCTION SIGN_ARRAY_3D_B
383#endif
384
385   !!======================================================================
386END MODULE lib_fortran
Note: See TracBrowser for help on using the repository browser.