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.
mppsumtam.F90 in branches/TAM_V3_2_2/NEMOTAM/OPATAM_SRC – NEMO

source: branches/TAM_V3_2_2/NEMOTAM/OPATAM_SRC/mppsumtam.F90 @ 3317

Last change on this file since 3317 was 2578, checked in by rblod, 13 years ago

first import of NEMOTAM 3.2.2
File size: 14.8 KB
Line 
1MODULE mppsumtam
2   !!======================================================================
3   !!                       ***  MODULE mppsumtam  ***
4   !! NEMOTAM: Summation of arrays across processors
5   !!======================================================================
6
7   !!----------------------------------------------------------------------
8   !!   mpp_sum_inter  : Interface for depending on nmppsum
9   !!   mpp_sum_simple : Simple sum
10   !!   loc_sum_inter  : Interface for depending on nmppsum (local version)
11   !!   loc_sum_simple : Simple sum (local version)
12   !!   loc_sum_indep  : Order independent MPP reproducible sum
13   !!----------------------------------------------------------------------
14   !! * Modules used   
15   USE par_kind, ONLY : &   ! Precision variables
16      & wp
17   USE dom_oce, ONLY : &    ! Ocean space and time domain variables
18      & nproc
19   USE par_oce, ONLY : &    ! Ocean parameters
20      & jpnij
21#if defined key_mpp_mpi
22   USE lib_mpp, ONLY : &    ! MPP library
23      & mpi_comm_opa,  &
24      & mpp_sum
25#endif
26   USE lib_mpp
27   USE mpp_tam
28   USE mppsum
29   USE in_out_manager
30
31   IMPLICIT NONE
32
33   !! * Routine accessibility
34   PRIVATE
35
36   PUBLIC &
37      & mpp_sum_inter,  & ! Interface for depending on nmppsum
38      & mpp_sum_simple, & ! Simple sum
39      & loc_sum_inter,  & ! Interface for depending on nmppsum
40      & loc_sum_simple, & ! Simple sum
41      & loc_sum_indep,  & ! Order independent local reproducible sum
42      & nmppsum           ! Summation control
43
44   PUBLIC mpp_sum_nfd
45
46   INTERFACE mpp_sum_nfd
47      MODULE PROCEDURE mpp_sum_nfd_3d,mpp_sum_nfd_4d 
48   END INTERFACE
49
50   !! * Module variables
51   INTEGER :: nmppsum = 2 ! Choice of MPP sum for
52                          ! (1 = simple)
53                          ! (2 = order independent)
54   
55CONTAINS
56
57   FUNCTION mpp_sum_inter( pval, kn )
58      !!----------------------------------------------------------------------
59      !!               ***  ROUTINE mpp_sum_inter ***
60      !!         
61      !! ** Purpose : Summation of arrays across processors
62      !!
63      !! ** Method  : Call either mpp_sum_simple or mpp_sum_indep depending on
64      !!              nmppsum
65      !!
66      !! ** Action  : This does only work for MPI.
67      !!              It does not work for SHMEM.
68      !!
69      !! References : http://www.mpi-forum.org
70      !!
71      !! History :
72      !!        !  07-07  (K. Mogensen)  Original code
73      !!----------------------------------------------------------------------
74      !! * Function return
75      REAL(wp) mpp_sum_inter
76      !! * Arguments
77      INTEGER, INTENT(IN) :: &
78         & kn
79      REAL(wp), DIMENSION(kn), INTENT(IN) :: &
80         & pval
81      !! * Local declarations
82     
83      IF ( nmppsum == 1 ) THEN
84         
85         mpp_sum_inter = mpp_sum_simple( pval, kn )
86         
87      ELSEIF ( nmppsum == 2 ) THEN
88
89         mpp_sum_inter = mpp_sum_indep( pval, kn )
90
91      ELSE
92
93         CALL ctl_stop( ' mpp_sum: Invalid nmppsum ' )
94
95      ENDIF
96
97   END FUNCTION mpp_sum_inter
98
99   FUNCTION mpp_sum_simple( pval, kn )
100      !!----------------------------------------------------------------------
101      !!               ***  ROUTINE mpp_sum_simple ***
102      !!         
103      !! ** Purpose : Summation of arrays across processors
104      !!
105      !! ** Method  : This routine is the naive implementation using
106      !!              mpi_allreduce for the sum of local data
107      !!
108      !! ** Action  : This does only work for MPI.
109      !!              It does not work for SHMEM.
110      !!
111      !! References : http://www.mpi-forum.org
112      !!
113      !! History :
114      !!        !  07-07  (K. Mogensen)  Original code
115      !!----------------------------------------------------------------------
116      !! * Function return
117      REAL(wp) mpp_sum_simple
118      !! * Arguments
119      INTEGER, INTENT(IN) :: &
120         & kn
121      REAL(wp), DIMENSION(kn), INTENT(IN) :: &
122         & pval
123      !! * Local declarations
124      REAL(wp) :: &
125         & zres
126      INTEGER :: &
127         & ierr
128      INTEGER :: &
129         & ji
130
131      ! Compute local sum   
132      zres = SUM( pval(1:kn) )   
133      ! Global sum
134
135      CALL mpp_sum( zres )
136
137      mpp_sum_simple = zres
138
139   END FUNCTION mpp_sum_simple
140
141   FUNCTION loc_sum_inter( pval, kn )
142      !!----------------------------------------------------------------------
143      !!               ***  ROUTINE loc_sum_inter ***
144      !!         
145      !! ** Purpose : Summation of arrays across processors
146      !!
147      !! ** Method  : Call either loc_sum_simple or loc_sum_indep depending on
148      !!              nmppsum
149      !!
150      !! ** Action  :
151      !!
152      !! References :
153      !!
154      !! History :
155      !!        !  07-09  (K. Mogensen)  Original code
156      !!----------------------------------------------------------------------
157      !! * Function return
158      REAL(wp) loc_sum_inter
159      !! * Arguments
160      INTEGER, INTENT(IN) :: &
161         & kn
162      REAL(wp), DIMENSION(kn), INTENT(IN) :: &
163         & pval
164      !! * Local declarations
165     
166      IF ( nmppsum == 1 ) THEN
167         
168         loc_sum_inter = loc_sum_simple( pval, kn )
169         
170      ELSEIF ( nmppsum == 2 ) THEN
171
172         loc_sum_inter = loc_sum_indep( pval, kn )
173
174      ELSE
175
176         CALL ctl_stop( ' loc_sum: Invalid nmppsum = ' )
177
178      ENDIF
179
180   END FUNCTION loc_sum_inter
181
182   FUNCTION loc_sum_simple( pval, kn )
183      !!----------------------------------------------------------------------
184      !!               ***  ROUTINE loc_sum_simple ***
185      !!         
186      !! ** Purpose : Summation of arrays across processors
187      !!
188      !! ** Method  : This routine is the naive implementation for the sum
189      !!              of local data
190      !!
191      !! ** Action  :
192      !!
193      !! References :
194      !!
195      !! History :
196      !!        !  07-09  (K. Mogensen)  Original code
197      !!----------------------------------------------------------------------
198      !! * Function return
199      REAL(wp) loc_sum_simple
200      !! * Arguments
201      INTEGER, INTENT(IN) :: &
202         & kn
203      REAL(wp), DIMENSION(kn), INTENT(IN) :: &
204         & pval
205      !! * Local declarations
206      REAL(wp) :: &
207         & zloc
208      INTEGER :: &
209         & ierr
210      INTEGER :: &
211         & ji
212
213      ! Compute local sum
214
215      zloc = SUM( pval(1:kn) )
216     
217      loc_sum_simple = zloc
218
219   END FUNCTION loc_sum_simple
220
221   FUNCTION loc_sum_indep( pval, kn )
222      !!----------------------------------------------------------------------
223      !!               ***  ROUTINE loc_sum_indep ***
224      !!         
225      !! ** Purpose : Sum all elements in the pval array in
226      !!              an accurate order-independent way.
227      !!
228      !! ** Method  : The code iterates the compensated summation until the
229      !!              result is guaranteed to be within 4*eps of the true sum.
230      !!              It then rounds the result to the nearest floating-point
231      !!              number whose last three bits are zero, thereby
232      !!              guaranteeing an order-independent result.
233      !!
234      !!              This version is used for local arrays. It is identical
235      !!              to the MPP version for clarity.
236      !!
237      !! ** Action  :
238      !!
239      !! References : M. Fisher (ECMWF): IFS code + personal communication
240      !!              The algorithm is based on Ogita et al. (2005)
241      !!              SIAM J. Sci. Computing, Vol.26, No.6, pp1955-1988.
242      !!              This is based in turn on an algorithm
243      !!              by Knuth (1969, seminumerical algorithms).
244      !!
245      !! History :
246      !!        !  07-09  (K. Mogensen)  Original code heavily based on IFS.
247      !!----------------------------------------------------------------------
248      !! * Function return
249      REAL(wp) loc_sum_indep
250      !! * Arguments
251      INTEGER, INTENT(IN) :: &
252         & kn
253      REAL(wp), DIMENSION(kn), INTENT(IN) :: &
254         & pval
255      !! * Local declarations
256      REAL(wp), DIMENSION(3) ::&
257         & zbuffl
258      REAL(wp), DIMENSION(:), ALLOCATABLE :: &
259         & zpsums, &
260         & zperrs, &
261         & zpcors, &
262         & zbuffg, &
263         & zp
264      REAL(wp) :: &
265         & zcorr,   &
266         & zerr,    &
267         & zolderr, &
268         & zbeta,   &
269         & zres
270      INTEGER :: &
271         & jj
272
273      ! Check that the the algorithm can work
274     
275      IF ( ( REAL( 2 * kn ) * EPSILON( zres ) ) >= 1.0 ) THEN
276
277         CALL ctl_stop( 'kn is too large to guarantee error bounds' )
278
279      ENDIF
280     
281      ALLOCATE( &
282         & zp(MAX(kn,1)),                & 
283         & zbuffg(1*SIZE(zbuffl)), &
284         & zpsums(1),              &
285         & zperrs(1),              &
286         & zpcors(1)               &
287         & )
288
289      zolderr = HUGE(zerr)
290
291      ! Copy the input array. This avoids some tricky indexing, at the
292      ! expense of some inefficency.
293
294      IF ( kn > 0 ) THEN
295
296         zp(:) = pval(:)
297
298      ELSE
299
300         zp(1) = 0.0_wp
301
302      ENDIF
303
304      k_loop: DO
305
306         ! Transform local arrays
307         
308         IF ( kn > 0 ) THEN
309
310            CALL comp_sum ( zp, kn, zcorr, zerr )
311
312         ENDIF
313
314         ! Gather partial sums and error bounds to all processors
315
316         zbuffl(1) = zp(MAX(kn,1))
317         
318         IF ( kn > 0 ) THEN
319
320            zbuffl(2) = zerr
321            zbuffl(3) = zcorr
322
323         ELSE
324           
325            zbuffl(2) = 0.0_wp
326            zbuffl(3) = 0.0_wp
327
328         ENDIF
329
330         zpsums(1) = zbuffl(1)
331         zperrs(1) = zbuffl(2)
332         zpcors(1) = zbuffl(3)
333
334         ! Transform partial sums
335
336         CALL comp_sum( zpsums, 1, zcorr, zerr )
337         zerr  = zerr  + SUM(zperrs)
338         zcorr = zcorr + SUM(zpcors)
339         
340         ! Calculate final result
341         
342         zres = zpsums(1) + zcorr
343
344         ! Calculate error bound. This is corollary 4.7 from Ogita et al.
345         ! (2005)
346         
347         zbeta = zerr *( REAL( 2*kn, wp ) * EPSILON(zres) ) &
348            &  /(1.0_wp - REAL( 2*kn, wp ) * EPSILON(zres) )
349
350         zerr = EPSILON(zres) * ABS(zres) &
351            & +(zbeta + ( 2.0_wp * EPSILON(zres) * EPSILON(zres) * ABS(zres) &
352            &            +3.0_wp * TINY(zres) ) )
353
354         ! Update the last element of the local array
355
356         zp(MAX(kn,1)) = zpsums(nproc+1)
357
358         ! Exit if the global error is small enough
359
360         IF ( zerr < 4.0_wp * SPACING(zres) ) EXIT k_loop
361         
362         ! Take appropriate action if ZRES cannot be sufficiently refined.
363         
364         IF (zerr >= zolderr) THEN
365
366            CALL ctl_stop( 'Failed to refine sum', &
367               &           'Warning: Possiblity of non-reproducible results' )
368
369         ENDIF
370
371         zolderr = zerr
372         
373      ENDDO k_loop
374
375      ! At this stage, we have guaranteed that ZRES less than 4*EPS
376      ! away from the exact sum. There are only four floating point
377      ! numbers in this range. So, if we find the nearest number that
378      ! has its last three bits zero, then we have a reproducible result.
379
380      loc_sum_indep = fround(zres)
381
382      DEALLOCATE( &
383         & zpcors, &
384         & zperrs, &
385         & zpsums, &
386         & zbuffg, &
387         & zp      &
388         & )
389
390   END FUNCTION loc_sum_indep
391
392#if defined key_mpp_mpi
393   FUNCTION mpp_sum_nfd_4d( pval, kn1, kn2, kn3, kn4, kcom )
394      !!----------------------------------------------------------------------
395      !!               ***  ROUTINE mpp_sum_nfd_4d ***
396      !!         
397      !! ** Purpose : Summation of arrays across processors (of kcom group)
398      !!
399      !! ** Method  : Pack 4d array to vector, sum vector and then unpack
400      !!
401      !! ** Action  : This does only work for MPI.
402      !!              It does not work for SHMEM.
403      !!
404      !! History :
405      !!        !  10-11  (F. Vigilant)  Original code
406      !!----------------------------------------------------------------------
407      !! * Function return
408      REAL(wp), DIMENSION(kn1, kn2, kn3, kn4) :: mpp_sum_nfd_4d
409      !! * Arguments
410      INTEGER, INTENT(IN) :: kn1, kn2, kn3, kn4
411      REAL(wp), DIMENSION(kn1, kn2, kn3, kn4), INTENT(IN) :: pval
412      INTEGER , INTENT( in ), OPTIONAL        :: kcom
413      !! * Local declarations
414      REAL(wp), DIMENSION(:), ALLOCATABLE     :: zvec
415      LOGICAL, DIMENSION(kn1, kn2, kn3, kn4)  :: zmask
416      REAL(wp), DIMENSION(kn1, kn2, kn3, kn4) :: zfd
417      INTEGER :: zdim
418     
419      zdim = kn1 * kn2 * kn3 * kn4
420      ALLOCATE( zvec(zdim) )
421     
422      zvec = PACK( pval(1:kn1,1:kn2,1:kn3,1:kn4),.TRUE.)
423
424      CALL mpp_sum( zvec, zdim, kcom )
425
426      zmask(:,:,:,:) = .TRUE.
427      mpp_sum_nfd_4d = UNPACK( zvec, zmask, zfd )
428     
429   END FUNCTION mpp_sum_nfd_4d
430
431   FUNCTION mpp_sum_nfd_3d( pval, kn1, kn2, kn3, kcom )
432      !!----------------------------------------------------------------------
433      !!               ***  ROUTINE mpp_sum_nfd_3d ***
434      !!         
435      !! ** Purpose : Summation of arrays across processors (of kcom group)
436      !!
437      !! ** Method  : Pack 3d array to vector, sum vector and then unpack
438      !!
439      !! ** Action  : This does only work for MPI.
440      !!              It does not work for SHMEM.
441      !!
442      !! History :
443      !!        !  10-11  (F. Vigilant)  Original code
444      !!----------------------------------------------------------------------
445      !! * Function return
446      REAL(wp), DIMENSION(kn1, kn2, kn3) :: mpp_sum_nfd_3d
447      !! * Arguments
448      INTEGER, INTENT(IN) ::kn1, kn2, kn3
449      REAL(wp), DIMENSION(kn1, kn2, kn3), INTENT(IN) :: pval
450      INTEGER , INTENT( in ), OPTIONAL        :: kcom
451      !! * Local declarations
452      REAL(wp), DIMENSION(:), ALLOCATABLE     :: zvec
453      LOGICAL, DIMENSION(kn1, kn2, kn3)  :: zmask
454      REAL(wp), DIMENSION(kn1, kn2, kn3) :: zfd
455      INTEGER :: zdim
456     
457      zdim = kn1 * kn2 * kn3
458      ALLOCATE( zvec(zdim) )
459     
460      zvec = PACK( pval(1:kn1,1:kn2,1:kn3),.TRUE.)
461
462      CALL mpp_sum( zvec, zdim, kcom )
463
464      zmask(:,:,:) = .TRUE.
465      mpp_sum_nfd_3d = UNPACK( zvec, zmask, zfd )
466     
467   END FUNCTION mpp_sum_nfd_3d
468#else
469   FUNCTION mpp_sum_nfd_4d( pval, kn1, kn2, kn3, kn4, kcom )
470      !! * Function return
471      REAL(wp), DIMENSION(kn1, kn2, kn3, kn4) :: mpp_sum_nfd_4d
472      !! * Arguments
473      INTEGER, INTENT(IN) :: kn1, kn2, kn3, kn4
474      REAL(wp), DIMENSION(kn1, kn2, kn3, kn4), INTENT(IN) :: pval
475      INTEGER , INTENT( in ), OPTIONAL        :: kcom
476      mpp_sum_nfd_4d = 0.0_wp
477      WRITE(*,*) 'mpp_sum_nfd_4d: You should not have seen this print! error?'
478   END FUNCTION mpp_sum_nfd_4d
479   FUNCTION mpp_sum_nfd_3d( pval, kn1, kn2, kn3, kcom )
480      !! * Function return
481      REAL(wp), DIMENSION(kn1, kn2, kn3) :: mpp_sum_nfd_3d
482      !! * Arguments
483      INTEGER, INTENT(IN) :: kn1, kn2, kn3
484      REAL(wp), DIMENSION(kn1, kn2, kn3), INTENT(IN) :: pval
485      INTEGER , INTENT( in ), OPTIONAL        :: kcom
486      mpp_sum_nfd_3d = 0.0_wp
487      WRITE(*,*) 'mpp_sum_nfd_3d: You should not have seen this print! error?'
488   END FUNCTION mpp_sum_nfd_3d
489#endif
490
491END MODULE mppsumtam
Note: See TracBrowser for help on using the repository browser.