FMS  2025.04
Flexible Modeling System
mpp_efp.F90
1 !***********************************************************************
2 !* Apache License 2.0
3 !*
4 !* This file is part of the GFDL Flexible Modeling System (FMS).
5 !*
6 !* Licensed under the Apache License, Version 2.0 (the "License");
7 !* you may not use this file except in compliance with the License.
8 !* You may obtain a copy of the License at
9 !*
10 !* http://www.apache.org/licenses/LICENSE-2.0
11 !*
12 !* FMS is distributed in the hope that it will be useful, but WITHOUT
13 !* WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied;
14 !* without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
15 !* PARTICULAR PURPOSE. See the License for the specific language
16 !* governing permissions and limitations under the License.
17 !***********************************************************************
18 !> @defgroup mpp_efp_mod mpp_efp_mod
19 !> @ingroup mpp
20 !> @brief This module provides interfaces to the non-domain-oriented communication
21 !! subroutines.
22 !!
23 !> Mainly includes interfaces and type definitions for reproducing operations with extended
24 !! fixed point data.
25 
26 !> @addtogroup mpp_efp_mod
27 !> @{
28 module mpp_efp_mod
29 
30 use mpp_mod, only : mpp_error, fatal, warning, note
31 use mpp_mod, only : mpp_pe, mpp_root_pe, mpp_npes
32 use mpp_mod, only : mpp_sum
33 use platform_mod
34 
35 implicit none ; private
36 
38 public :: mpp_efp_plus, mpp_efp_minus, mpp_efp_to_real, mpp_real_to_efp, mpp_efp_real_diff
39 public :: operator(+), operator(-), assignment(=)
40 public :: mpp_query_efp_overflow_error, mpp_reset_efp_overflow_error
41 
42 integer, parameter :: NUMBIT = 46 !< number of bits used in the 64-bit signed integer representation.
43 integer, parameter :: numint = 6 !< The number of long integers to use to represent
44  !! a real number.
45 
46 integer(i8_kind), parameter :: prec=2_8**numbit !< The precision of each integer.
47 real(r8_kind), parameter :: r_prec=2.0_8**numbit !< A real version of prec.
48 real(r8_kind), parameter :: i_prec=1.0_8/(2.0_8**numbit) !< The inverse of prec.
49 integer, parameter :: max_count_prec=2**(63-numbit)-1 !< The number of values that can be added together
50  !! with the current value of prec before there will
51  !! be roundoff problems.
52 
53 real(r8_kind), parameter, dimension(NUMINT) :: &
54  pr = (/ r_prec**2, r_prec, 1.0_8, 1.0_8/r_prec, 1.0_8/r_prec**2, 1.0_8/r_prec**3 /)
55 real(r8_kind), parameter, dimension(NUMINT) :: &
56  i_pr = (/ 1.0_8/r_prec**2, 1.0_8/r_prec, 1.0_8, r_prec, r_prec**2, r_prec**3 /)
57 
58 logical :: overflow_error = .false., nan_error = .false.
59 logical :: debug = .false. !< Making this true enables debugging output.
60 
61 !> @}
62 
63 !> This interface uses a conversion to an integer representation
64 !! of real numbers to give order-invariant sums that will reproduce
65 !! across PE count.
66 !!
67 !! This idea comes from R. Hallberg and A. Adcroft.
68 !!
69 !> @ingroup mpp_efp_mod
71  module procedure mpp_reproducing_sum_r8_2d
72  module procedure mpp_reproducing_sum_r8_3d
73  module procedure mpp_reproducing_sum_r4_2d
74 end interface mpp_reproducing_sum
75 
76 !> The Extended Fixed Point (mpp_efp) type provides a public interface for doing
77 !! sums and taking differences with this type.
78 !> @ingroup mpp_efp_mod
79 type, public :: mpp_efp_type
80  private
81  integer(i8_kind), dimension(NUMINT) :: v
82 end type mpp_efp_type
83 
84 
85 !> Operator override interface for mpp_efp_type
86 !> @ingroup mpp_efp_mod
87 interface operator (+); module procedure mpp_efp_plus ; end interface
88 !> Operator override interface for mpp_efp_type
89 !> @ingroup mpp_efp_mod
90 interface operator (-); module procedure mpp_efp_minus ; end interface
91 !> Assignment override interface for mpp_efp_type
92 !> @ingroup mpp_efp_mod
93 interface assignment(=); module procedure mpp_efp_assign ; end interface
94 
95 !> @addtogroup mpp_efp_mod
96 !> @{
97 
98 contains
99 
100 !> @brief Calculates a reproducing sum for a 2D, 8-byte real array
101 function mpp_reproducing_sum_r8_2d(array, isr, ier, jsr, jer, EFP_sum, reproducing, &
102  overflow_check, err) result(sum)
103  real(r8_kind), dimension(:,:), intent(in) :: array
104  integer, optional, intent(in) :: isr, ier, jsr, jer
105  type(mpp_efp_type), optional, intent(out) :: efp_sum
106  logical, optional, intent(in) :: reproducing
107  logical, optional, intent(in) :: overflow_check
108  integer, optional, intent(out) :: err
109  real(r8_kind) :: sum !< Result
110 
111  integer(i8_kind), dimension(NUMINT) :: ints_sum
112  integer(i8_kind) :: ival, prec_error
113  real(r8_kind) :: rsum(1), rs
114  real(r8_kind) :: max_mag_term
115  logical :: repro, over_check
116  character(len=256) :: mesg
117  integer :: i, j, n, is, ie, js, je, sgn
118 
119  if (mpp_npes() > max_count_prec) call mpp_error(fatal, &
120  "mpp_reproducing_sum: Too many processors are being used for the value of "//&
121  "prec. Reduce prec to (2^63-1)/mpp_npes.")
122 
123  prec_error = (2_8**62 + (2_8**62 - 1)) / mpp_npes()
124 
125  is = 1 ; ie = size(array,1) ; js = 1 ; je = size(array,2 )
126  if (present(isr)) then
127  if (isr < is) call mpp_error(fatal, &
128  "Value of isr too small in mpp_reproducing_sum_2d.")
129  is = isr
130  endif
131  if (present(ier)) then
132  if (ier > ie) call mpp_error(fatal, &
133  "Value of ier too large in mpp_reproducing_sum_2d.")
134  ie = ier
135  endif
136  if (present(jsr)) then
137  if (jsr < js) call mpp_error(fatal, &
138  "Value of jsr too small in mpp_reproducing_sum_2d.")
139  js = jsr
140  endif
141  if (present(jer)) then
142  if (jer > je) call mpp_error(fatal, &
143  "Value of jer too large in mpp_reproducing_sum_2d.")
144  je = jer
145  endif
146 
147  repro = .true. ; if (present(reproducing)) repro = reproducing
148  over_check = .true. ; if (present(overflow_check)) over_check = overflow_check
149 
150  if (repro) then
151  overflow_error = .false. ; nan_error = .false. ; max_mag_term = 0.0
152  ints_sum(:) = 0
153  if (over_check) then
154  if ((je+1-js)*(ie+1-is) < max_count_prec) then
155  do j=js,je ; do i=is,ie
156  call increment_ints_faster(ints_sum, array(i,j), max_mag_term);
157  enddo ; enddo
158  call carry_overflow(ints_sum, prec_error)
159  elseif ((ie+1-is) < max_count_prec) then
160  do j=js,je
161  do i=is,ie
162  call increment_ints_faster(ints_sum, array(i,j), max_mag_term);
163  enddo
164  call carry_overflow(ints_sum, prec_error)
165  enddo
166  else
167  do j=js,je ; do i=is,ie
168  call increment_ints(ints_sum, real_to_ints(array(i,j), prec_error), &
169  prec_error);
170  enddo ; enddo
171  endif
172  else
173  do j=js,je ; do i=is,ie
174  sgn = 1 ; if (array(i,j)<0.0) sgn = -1
175  rs = abs(array(i,j))
176  do n=1,numint
177  ival = int(rs*i_pr(n), 8)
178  rs = rs - ival*pr(n)
179  ints_sum(n) = ints_sum(n) + sgn*ival
180  enddo
181  enddo ; enddo
182  call carry_overflow(ints_sum, prec_error)
183  endif
184 
185  if (present(err)) then
186  err = 0
187  if (overflow_error) &
188  err = err+2
189  if (nan_error) &
190  err = err+4
191  if (err > 0) then ; do n=1,numint ; ints_sum(n) = 0 ; enddo ; endif
192  else
193  if (nan_error) then
194  call mpp_error(fatal, "NaN in input field of mpp_reproducing_sum(_2d), this indicates numerical instability")
195  endif
196  if (abs(max_mag_term) >= prec_error*pr(1)) then
197  write(mesg, '(ES13.5)') max_mag_term
198  call mpp_error(fatal,"Overflow in mpp_reproducing_sum(_2d) conversion of "//trim(mesg))
199  endif
200  if (overflow_error) then
201  call mpp_error(fatal, "Overflow in mpp_reproducing_sum(_2d).")
202  endif
203  endif
204 
205  call mpp_sum(ints_sum, numint)
206 
207  call regularize_ints(ints_sum)
208  sum = ints_to_real(ints_sum)
209  else
210  rsum(1) = 0.0
211  do j=js,je ; do i=is,ie
212  rsum(1) = rsum(1) + array(i,j);
213  enddo ; enddo
214  call mpp_sum(rsum,1)
215  sum = rsum(1)
216 
217  if (present(err)) then ; err = 0 ; endif
218 
219  if (debug .or. present(efp_sum)) then
220  overflow_error = .false.
221  ints_sum = real_to_ints(sum, prec_error, overflow_error)
222  if (overflow_error) then
223  if (present(err)) then
224  err = err + 2
225  else
226  write(mesg, '(ES13.5)') sum
227  call mpp_error(fatal,"Repro_sum_2d: Overflow in real_to_ints conversion of "//trim(mesg))
228  endif
229  endif
230  endif
231  endif
232 
233  if (present(efp_sum)) efp_sum%v(:) = ints_sum(:)
234 
235  if (debug) then
236  write(mesg,'("2d RS: ", ES24.16, 6 Z17.16)') sum, ints_sum(1:numint)
237  if(mpp_pe() == mpp_root_pe()) call mpp_error(note, mesg)
238  endif
239 
240 end function mpp_reproducing_sum_r8_2d
241 
242 function mpp_reproducing_sum_r4_2d(array, isr, ier, jsr, jer, EFP_sum, reproducing, &
243  overflow_check, err) result(sum)
244  real(r4_kind), dimension(:,:), intent(in) :: array
245  integer, optional, intent(in) :: isr, ier, jsr, jer
246  type(mpp_efp_type), optional, intent(out) :: efp_sum
247  logical, optional, intent(in) :: reproducing
248  logical, optional, intent(in) :: overflow_check
249  integer, optional, intent(out) :: err
250  real(r4_kind) :: sum !< Result
251 
252  real(r8_kind) :: array_r8(size(array,1), size(array,2))
253 
254  array_r8 = array
255 
256  sum = real(mpp_reproducing_sum_r8_2d(array_r8, isr, ier, jsr, jer, efp_sum, reproducing, &
257  overflow_check, err), r4_kind)
258 
259  return
260 
261 end function mpp_reproducing_sum_r4_2d
262 
263 !> @brief Reproducing sum for 3d arrays of 8-bit reals
264 !!
265 !> This function uses a conversion to an integer representation
266 !! of real numbers to give order-invariant sums that will reproduce
267 !! across PE count. This idea comes from R. Hallberg and A. Adcroft.
268 function mpp_reproducing_sum_r8_3d(array, isr, ier, jsr, jer, sums, EFP_sum, err) &
269  result(sum)
270  real(r8_kind), dimension(:,:,:), intent(in) :: array
271  integer, optional, intent(in) :: isr, ier, jsr, jer
272  real(r8_kind), dimension(:), optional, intent(out) :: sums
273  type(mpp_efp_type), optional, intent(out) :: efp_sum
274  integer, optional, intent(out) :: err
275  real(r8_kind) :: sum !< Result
276 
277  real(r8_kind) :: max_mag_term
278  integer(i8_kind), dimension(NUMINT) :: ints_sum
279  integer(i8_kind), dimension(NUMINT,size(array,3)) :: ints_sums
280  integer(i8_kind) :: prec_error
281  character(len=256) :: mesg
282  integer :: i, j, k, is, ie, js, je, ke, isz, jsz, n
283 
284  if (mpp_npes() > max_count_prec) call mpp_error(fatal, &
285  "mpp_reproducing_sum: Too many processors are being used for the value of "//&
286  "prec. Reduce prec to (2^63-1)/mpp_npes.")
287 
288  prec_error = (2_8**62 + (2_8**62 - 1)) / mpp_npes()
289  max_mag_term = 0.0
290 
291  is = 1 ; ie = size(array,1) ; js = 1 ; je = size(array,2) ; ke = size(array,3)
292  if (present(isr)) then
293  if (isr < is) call mpp_error(fatal, &
294  "Value of isr too small in mpp_reproducing_sum(_3d).")
295  is = isr
296  endif
297  if (present(ier)) then
298  if (ier > ie) call mpp_error(fatal, &
299  "Value of ier too large in mpp_reproducing_sum(_3d).")
300  ie = ier
301  endif
302  if (present(jsr)) then
303  if (jsr < js) call mpp_error(fatal, &
304  "Value of jsr too small in mpp_reproducing_sum(_3d).")
305  js = jsr
306  endif
307  if (present(jer)) then
308  if (jer > je) call mpp_error(fatal, &
309  "Value of jer too large in mpp_reproducing_sum(_3d).")
310  je = jer
311  endif
312  jsz = je+1-js; isz = ie+1-is
313 
314  if (present(sums)) then
315  if (size(sums) > ke) call mpp_error(fatal, "Sums is smaller than "//&
316  "the vertical extent of array in mpp_reproducing_sum(_3d).")
317  ints_sums(:,:) = 0
318  overflow_error = .false. ; nan_error = .false. ; max_mag_term = 0.0
319  if (jsz*isz < max_count_prec) then
320  do k=1,ke
321  do j=js,je ; do i=is,ie
322  call increment_ints_faster(ints_sums(:,k), array(i,j,k), max_mag_term);
323  enddo ; enddo
324  call carry_overflow(ints_sums(:,k), prec_error)
325  enddo
326  elseif (isz < max_count_prec) then
327  do k=1,ke ; do j=js,je
328  do i=is,ie
329  call increment_ints_faster(ints_sums(:,k), array(i,j,k), max_mag_term);
330  enddo
331  call carry_overflow(ints_sums(:,k), prec_error)
332  enddo ; enddo
333  else
334  do k=1,ke ; do j=js,je ; do i=is,ie
335  call increment_ints(ints_sums(:,k), &
336  real_to_ints(array(i,j,k), prec_error), prec_error);
337  enddo ; enddo ; enddo
338  endif
339  if (present(err)) then
340  err = 0
341  if (abs(max_mag_term) >= prec_error*pr(1)) err = err+1
342  if (overflow_error) err = err+2
343  if (nan_error) err = err+2
344  if (err > 0) then ; do k=1,ke ; do n=1,numint ; ints_sums(n,k) = 0 ; enddo ; enddo ; endif
345  else
346  if (nan_error) call mpp_error(fatal, &
347  "NaN in input field of mpp_reproducing_sum(_3d), this indicates numerical instability")
348  if (abs(max_mag_term) >= prec_error*pr(1)) then
349  write(mesg, '(ES13.5)') max_mag_term
350  call mpp_error(fatal,"Overflow in mpp_reproducing_sum(_3d) conversion of "//trim(mesg))
351  endif
352  if (overflow_error) call mpp_error(fatal, "Overflow in mpp_reproducing_sum(_3d).")
353  endif
354 
355  call mpp_sum(ints_sums(:,1:ke), numint*ke)
356 
357  sum = 0.0
358  do k=1,ke
359  call regularize_ints(ints_sums(:,k))
360  sums(k) = ints_to_real(ints_sums(:,k))
361  sum = sum + sums(k)
362  enddo
363 
364  if (present(efp_sum)) then
365  efp_sum%v(:) = 0
366  do k=1,ke ; call increment_ints(efp_sum%v(:), ints_sums(:,k)) ; enddo
367  endif
368 
369  if (debug) then
370  do n=1,numint ; ints_sum(n) = 0 ; enddo
371  do k=1,ke ; do n=1,numint ; ints_sum(n) = ints_sum(n) + ints_sums(n,k) ; enddo ; enddo
372  write(mesg,'("3D RS: ", ES24.16, 6 Z17.16)') sum, ints_sum(1:numint)
373  if(mpp_pe()==mpp_root_pe()) call mpp_error(note, mesg)
374  endif
375  else
376  ints_sum(:) = 0
377  overflow_error = .false. ; nan_error = .false. ; max_mag_term = 0.0
378  if (jsz*isz < max_count_prec) then
379  do k=1,ke
380  do j=js,je ; do i=is,ie
381  call increment_ints_faster(ints_sum, array(i,j,k), max_mag_term);
382  enddo ; enddo
383  call carry_overflow(ints_sum, prec_error)
384  enddo
385  elseif (isz < max_count_prec) then
386  do k=1,ke ; do j=js,je
387  do i=is,ie
388  call increment_ints_faster(ints_sum, array(i,j,k), max_mag_term);
389  enddo
390  call carry_overflow(ints_sum, prec_error)
391  enddo ; enddo
392  else
393  do k=1,ke ; do j=js,je ; do i=is,ie
394  call increment_ints(ints_sum, real_to_ints(array(i,j,k), prec_error), &
395  prec_error);
396  enddo ; enddo ; enddo
397  endif
398  if (present(err)) then
399  err = 0
400  if (abs(max_mag_term) >= prec_error*pr(1)) err = err+1
401  if (overflow_error) err = err+2
402  if (nan_error) err = err+2
403  if (err > 0) then ; do n=1,numint ; ints_sum(n) = 0 ; enddo ; endif
404  else
405  if (nan_error) call mpp_error(fatal, &
406  "NaN in input field of mpp_reproducing_sum(_3d), this indicates numerical instability")
407  if (abs(max_mag_term) >= prec_error*pr(1)) then
408  write(mesg, '(ES13.5)') max_mag_term
409  call mpp_error(fatal,"Overflow in mpp_reproducing_sum(_3d) conversion of "//trim(mesg))
410  endif
411  if (overflow_error) call mpp_error(fatal, "Overflow in mpp_reproducing_sum(_3d).")
412  endif
413 
414  call mpp_sum(ints_sum, numint)
415 
416  call regularize_ints(ints_sum)
417  sum = ints_to_real(ints_sum)
418 
419  if (present(efp_sum)) efp_sum%v(:) = ints_sum(:)
420 
421  if (debug) then
422  write(mesg,'("3d RS: ", ES24.16, 6 Z17.16)') sum, ints_sum(1:numint)
423  if(mpp_pe()==mpp_root_pe()) call mpp_error(note, mesg)
424  endif
425  endif
426 
427 end function mpp_reproducing_sum_r8_3d
428 
429 !> @brief This function converts a real number to an equivalent representation
430 !! using several long integers.
431 function real_to_ints(r, prec_error, overflow) result(ints)
432  real(r8_kind), intent(in) :: r
433  integer(i8_kind), optional, intent(in) :: prec_error
434  logical, optional, intent(inout) :: overflow
435  integer(i8_kind), dimension(NUMINT) :: ints
436 
437  real(r8_kind) :: rs
438  character(len=80) :: mesg
439  integer(i8_kind) :: ival, prec_err
440  integer :: sgn, i
441 
442  prec_err = prec ; if (present(prec_error)) prec_err = prec_error
443  ints(:) = 0_8
444  if ((r >= 1e30) .eqv. (r < 1e30)) then ; nan_error = .true. ; return ; endif
445 
446  sgn = 1 ; if (r<0.0) sgn = -1
447  rs = abs(r)
448 
449  if (present(overflow)) then
450  if (.not.(rs < prec_err*pr(1))) overflow = .true.
451  if ((r >= 1e30) .eqv. (r < 1e30)) overflow = .true.
452  elseif (.not.(rs < prec_err*pr(1))) then
453  write(mesg, '(ES13.5)') r
454  call mpp_error(fatal,"Overflow in real_to_ints conversion of "//trim(mesg))
455  endif
456 
457  do i=1,numint
458  ival = int(rs*i_pr(i), 8)
459  rs = rs - ival*pr(i)
460  ints(i) = sgn*ival
461  enddo
462 
463 end function real_to_ints
464 
465 !> @brief This function reverses the conversion in real_to_ints.
466 function ints_to_real(ints) result(r)
467  integer(i8_kind), dimension(NUMINT), intent(in) :: ints
468  real(r8_kind) :: r
469 
470  integer :: i
471 
472  r = 0.0
473  do i=1,numint ; r = r + pr(i)*ints(i) ; enddo
474 end function ints_to_real
475 
476 !> @brief This subroutine increments a number with another, both using the integer
477 !! representation in real_to_ints.
478 subroutine increment_ints(int_sum, int2, prec_error)
479  integer(i8_kind), dimension(NUMINT), intent(inout) :: int_sum
480  integer(i8_kind), dimension(NUMINT), intent(in) :: int2
481  integer(i8_kind), optional, intent(in) :: prec_error
482 
483  integer :: i
484 
485  do i=numint,2,-1
486  int_sum(i) = int_sum(i) + int2(i)
487  ! Carry the local overflow.
488  if (int_sum(i) > prec) then
489  int_sum(i) = int_sum(i) - prec
490  int_sum(i-1) = int_sum(i-1) + 1
491  elseif (int_sum(i) < -prec) then
492  int_sum(i) = int_sum(i) + prec
493  int_sum(i-1) = int_sum(i-1) - 1
494  endif
495  enddo
496  int_sum(1) = int_sum(1) + int2(1)
497  if (present(prec_error)) then
498  if (abs(int_sum(1)) > prec_error) overflow_error = .true.
499  else
500  if (abs(int_sum(1)) > prec) overflow_error = .true.
501  endif
502 
503 end subroutine increment_ints
504 
505 !> @brief This subroutine increments a number with another, both using the integer
506 !! representation in real_to_ints, but without doing any carrying of overflow.
507 !! The entire operation is embedded in a single call for greater speed.
508 subroutine increment_ints_faster(int_sum, r, max_mag_term)
509  integer(i8_kind), dimension(NUMINT), intent(inout) :: int_sum
510  real(r8_kind), intent(in) :: r
511  real(r8_kind), intent(inout) :: max_mag_term
512 
513  real(r8_kind) :: rs
514  integer(i8_kind) :: ival
515  integer :: sgn, i
516 
517  if ((r >= 1e30) .eqv. (r < 1e30)) then ; nan_error = .true. ; return ; endif
518  sgn = 1 ; if (r<0.0) sgn = -1
519  rs = abs(r)
520  if (rs > abs(max_mag_term)) max_mag_term = r
521 
522  do i=1,numint
523  ival = int(rs*i_pr(i), 8)
524  rs = rs - ival*pr(i)
525  int_sum(i) = int_sum(i) + sgn*ival
526  enddo
527 
528 end subroutine increment_ints_faster
529 
530 !> @brief This subroutine handles carrying of the overflow.
531 subroutine carry_overflow(int_sum, prec_error)
532  integer(i8_kind), dimension(NUMINT), intent(inout) :: int_sum
533  integer(i8_kind), intent(in) :: prec_error
534 
535  integer :: i, num_carry
536 
537  do i=numint,2,-1 ; if (abs(int_sum(i)) > prec) then
538  num_carry = int(int_sum(i) * i_prec)
539  int_sum(i) = int_sum(i) - num_carry*prec
540  int_sum(i-1) = int_sum(i-1) + num_carry
541  endif ; enddo
542  if (abs(int_sum(1)) > prec_error) then
543  overflow_error = .true.
544  endif
545 
546 end subroutine carry_overflow
547 
548 !> @brief This subroutine carries the overflow, and then makes sure that
549 !! all integers are of the same sign as the overall value.
550 subroutine regularize_ints(int_sum)
551  integer(i8_kind), dimension(NUMINT), intent(inout) :: int_sum
552 
553  logical :: positive
554  integer :: i, num_carry
555 
556  do i=numint,2,-1 ; if (abs(int_sum(i)) > prec) then
557  num_carry = int(int_sum(i) * i_prec)
558  int_sum(i) = int_sum(i) - num_carry*prec
559  int_sum(i-1) = int_sum(i-1) + num_carry
560  endif ; enddo
561 
562  ! Determine the sign of the final number.
563  positive = .true.
564  do i=1,numint
565  if (abs(int_sum(i)) > 0) then
566  if (int_sum(i) < 0) positive = .false.
567  exit
568  endif
569  enddo
570 
571  if (positive) then
572  do i=numint,2,-1 ; if (int_sum(i) < 0) then
573  int_sum(i) = int_sum(i) + prec
574  int_sum(i-1) = int_sum(i-1) - 1
575  endif ; enddo
576  else
577  do i=numint,2,-1 ; if (int_sum(i) > 0) then
578  int_sum(i) = int_sum(i) - prec
579  int_sum(i-1) = int_sum(i-1) + 1
580  endif ; enddo
581  endif
582 
583 end subroutine regularize_ints
584 
585 function mpp_query_efp_overflow_error()
586  logical :: mpp_query_efp_overflow_error
587  mpp_query_efp_overflow_error = overflow_error
588 end function mpp_query_efp_overflow_error
589 
590 subroutine mpp_reset_efp_overflow_error()
591  overflow_error = .false.
592 end subroutine mpp_reset_efp_overflow_error
593 
594 function mpp_efp_plus(EFP1, EFP2)
595  type(mpp_efp_type) :: mpp_efp_plus
596  type(mpp_efp_type), intent(in) :: EFP1, EFP2
597 
598  mpp_efp_plus = efp1
599 
600  call increment_ints(mpp_efp_plus%v(:), efp2%v(:))
601 end function mpp_efp_plus
602 
603 function mpp_efp_minus(EFP1, EFP2)
604  type(mpp_efp_type) :: mpp_efp_minus
605  type(mpp_efp_type), intent(in) :: EFP1, EFP2
606  integer :: i
607 
608  do i=1,numint ; mpp_efp_minus%v(i) = -1*efp2%v(i) ; enddo
609 
610  call increment_ints(mpp_efp_minus%v(:), efp1%v(:))
611 end function mpp_efp_minus
612 
613 !> @brief This subroutine assigns all components of the extended fixed point type
614 !! variable on the RHS (EFP2) to the components of the variable on the LHS
615 !! (EFP1).
616 subroutine mpp_efp_assign(EFP1, EFP2)
617  type(mpp_efp_type), intent(out) :: EFP1
618  type(mpp_efp_type), intent(in) :: EFP2
619  integer i
620 
621  do i=1,numint ; efp1%v(i) = efp2%v(i) ; enddo
622 end subroutine mpp_efp_assign
623 
624 function mpp_efp_to_real(EFP1)
625  type(mpp_efp_type), intent(inout) :: EFP1
626  real(r8_kind) :: mpp_efp_to_real
627 
628  call regularize_ints(efp1%v)
629  mpp_efp_to_real = ints_to_real(efp1%v)
630 end function mpp_efp_to_real
631 
632 function mpp_efp_real_diff(EFP1, EFP2)
633  type(mpp_efp_type), intent(in) :: EFP1, EFP2
634  real(r8_kind) :: mpp_efp_real_diff
635 
636  type(mpp_efp_type) :: EFP_diff
637 
638  efp_diff = efp1 - efp2
639  mpp_efp_real_diff = mpp_efp_to_real(efp_diff)
640 
641 end function mpp_efp_real_diff
642 
643 function mpp_real_to_efp(val, overflow)
644  real(r8_kind), intent(in) :: val
645  logical, optional, intent(inout) :: overflow
646  type(mpp_efp_type) :: mpp_real_to_efp
647 
648  logical :: over
649  character(len=80) :: mesg
650 
651  if (present(overflow)) then
652  mpp_real_to_efp%v(:) = real_to_ints(val, overflow=overflow)
653  else
654  over = .false.
655  mpp_real_to_efp%v(:) = real_to_ints(val, overflow=over)
656  if (over) then
657  write(mesg, '(ES13.5)') val
658  call mpp_error(fatal,"Overflow in mpp_real_to_efp conversion of "//trim(mesg))
659  endif
660  endif
661 
662 end function mpp_real_to_efp
663 
664 !> This subroutine does a sum across PEs of a list of EFP variables,
665 !! returning the sums in place, with all overflows carried.
666 subroutine mpp_efp_list_sum_across_pes(EFPs, nval, errors)
667  type(mpp_efp_type), dimension(:), intent(inout) :: efps
668  integer, intent(in) :: nval
669  logical, dimension(:), optional, intent(out) :: errors
670 
671  integer(i8_kind), dimension(NUMINT,nval) :: ints
672  integer(i8_kind) :: prec_error
673  logical :: error_found
674  character(len=256) :: mesg
675  integer :: i, n
676 
677  if (mpp_npes() > max_count_prec) call mpp_error(fatal, &
678  "mpp_efp_list_sum_across_PEs: Too many processors are being used for the value of "//&
679  "prec. Reduce prec to (2^63-1)/mpp_npes.")
680 
681  prec_error = (2_8**62 + (2_8**62 - 1)) / mpp_npes()
682  ! overflow_error is an overflow error flag for the whole module.
683  overflow_error = .false. ; error_found = .false.
684 
685  do i=1,nval ; do n=1,numint ; ints(n,i) = efps(i)%v(n) ; enddo ; enddo
686 
687  call mpp_sum(ints(:,:), numint*nval)
688 
689  if (present(errors)) errors(:) = .false.
690  do i=1,nval
691  overflow_error = .false.
692  call carry_overflow(ints(:,i), prec_error)
693  do n=1,numint ; efps(i)%v(n) = ints(n,i) ; enddo
694  if (present(errors)) errors(i) = overflow_error
695  if (overflow_error) then
696  write (mesg,'("mpp_efp_list_sum_across_PEs error at ",i6," val was ",ES13.6, ", prec_error = ",ES13.6)') &
697  i, mpp_efp_to_real(efps(i)), real(prec_error)
698  if(mpp_pe()==mpp_root_pe()) call mpp_error(warning, mesg)
699  endif
700  error_found = error_found .or. overflow_error
701  enddo
702  if (error_found .and. .not.(present(errors))) then
703  call mpp_error(fatal, "Overflow in mpp_efp_list_sum_across_PEs.")
704  endif
705 
706 end subroutine mpp_efp_list_sum_across_pes
707 
708 end module mpp_efp_mod
709 !> @}
710 ! close documentation grouping
integer, parameter max_count_prec
The number of values that can be added together with the current value of prec before there will be r...
Definition: mpp_efp.F90:49
subroutine increment_ints_faster(int_sum, r, max_mag_term)
This subroutine increments a number with another, both using the integer representation in real_to_in...
Definition: mpp_efp.F90:509
logical debug
Making this true enables debugging output.
Definition: mpp_efp.F90:59
subroutine carry_overflow(int_sum, prec_error)
This subroutine handles carrying of the overflow.
Definition: mpp_efp.F90:532
real(r8_kind), parameter i_prec
The inverse of prec.
Definition: mpp_efp.F90:48
integer(i8_kind) function, dimension(numint) real_to_ints(r, prec_error, overflow)
This function converts a real number to an equivalent representation using several long integers.
Definition: mpp_efp.F90:432
real(r8_kind) function mpp_reproducing_sum_r8_2d(array, isr, ier, jsr, jer, EFP_sum, reproducing, overflow_check, err)
Calculates a reproducing sum for a 2D, 8-byte real array.
Definition: mpp_efp.F90:103
real(r8_kind), parameter r_prec
A real version of prec.
Definition: mpp_efp.F90:47
integer, parameter numint
The number of long integers to use to represent a real number.
Definition: mpp_efp.F90:43
integer(i8_kind), parameter prec
The precision of each integer.
Definition: mpp_efp.F90:46
real(r8_kind) function mpp_reproducing_sum_r8_3d(array, isr, ier, jsr, jer, sums, EFP_sum, err)
Reproducing sum for 3d arrays of 8-bit reals.
Definition: mpp_efp.F90:270
subroutine increment_ints(int_sum, int2, prec_error)
This subroutine increments a number with another, both using the integer representation in real_to_in...
Definition: mpp_efp.F90:479
real(r4_kind) function mpp_reproducing_sum_r4_2d(array, isr, ier, jsr, jer, EFP_sum, reproducing, overflow_check, err)
Definition: mpp_efp.F90:244
real(r8_kind) function ints_to_real(ints)
This function reverses the conversion in real_to_ints.
Definition: mpp_efp.F90:467
subroutine regularize_ints(int_sum)
This subroutine carries the overflow, and then makes sure that all integers are of the same sign as t...
Definition: mpp_efp.F90:551
subroutine, public mpp_efp_list_sum_across_pes(EFPs, nval, errors)
This subroutine does a sum across PEs of a list of EFP variables, returning the sums in place,...
Definition: mpp_efp.F90:667
This interface uses a conversion to an integer representation of real numbers to give order-invariant...
Definition: mpp_efp.F90:70
The Extended Fixed Point (mpp_efp) type provides a public interface for doing sums and taking differe...
Definition: mpp_efp.F90:79
integer function mpp_npes()
Returns processor count for current pelist.
Definition: mpp_util.inc:420
integer function mpp_pe()
Returns processor ID.
Definition: mpp_util.inc:406
Error handler.
Definition: mpp.F90:381
Reduction operation.
Definition: mpp.F90:596