30 use mpp_mod, 
only : 
mpp_error, fatal, warning, note
 
   35 implicit none ; 
private 
   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
 
   42 integer, 
parameter :: NUMBIT = 46
 
   46 integer(i8_kind), 
parameter :: 
prec=2_8**numbit 
 
   47 real(r8_kind), 
parameter :: 
r_prec=2.0_8**numbit 
 
   48 real(r8_kind), 
parameter :: 
i_prec=1.0_8/(2.0_8**numbit) 
 
   53 real(r8_kind), 
parameter, 
dimension(NUMINT) :: &
 
   55 real(r8_kind), 
parameter, 
dimension(NUMINT) :: &
 
   58 logical :: overflow_error = .false., nan_error = .false.
 
   81   integer(i8_kind), 
dimension(NUMINT) :: v
 
   87 interface operator (+); module 
procedure mpp_efp_plus  ; end interface
 
   90 interface operator (-); module 
procedure mpp_efp_minus ; end interface
 
   93 interface assignment(=); module 
procedure mpp_efp_assign ; end interface
 
  102                             overflow_check, err) 
result(sum)
 
  103   real(r8_kind), 
dimension(:,:), 
intent(in) :: array
 
  104   integer,        
optional,          
intent(in) :: isr, ier, jsr, jer
 
  106   logical,        
optional,          
intent(in) :: reproducing
 
  107   logical,        
optional,          
intent(in) :: overflow_check
 
  108   integer,        
optional,         
intent(out) :: err
 
  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
 
  120     "mpp_reproducing_sum: Too many processors are being used for the value of "//&
 
  121     "prec.  Reduce prec to (2^63-1)/mpp_npes.")
 
  123   prec_error = (2_8**62 + (2_8**62 - 1)) / 
mpp_npes()
 
  125   is = 1 ; ie = 
size(array,1) ; js = 1 ; je = 
size(array,2 )
 
  126   if (
present(isr)) 
then 
  128       "Value of isr too small in mpp_reproducing_sum_2d.")
 
  131   if (
present(ier)) 
then 
  133       "Value of ier too large in mpp_reproducing_sum_2d.")
 
  136   if (
present(jsr)) 
then 
  138       "Value of jsr too small in mpp_reproducing_sum_2d.")
 
  141   if (
present(jer)) 
then 
  143       "Value of jer too large in mpp_reproducing_sum_2d.")
 
  147   repro = .true. ; 
if (
present(reproducing)) repro = reproducing
 
  148   over_check = .true. ; 
if (
present(overflow_check)) over_check = overflow_check
 
  151     overflow_error = .false. ; nan_error = .false. ; max_mag_term = 0.0
 
  155         do j=js,je ; 
do i=is,ie
 
  167         do j=js,je ; 
do i=is,ie
 
  173       do j=js,je ; 
do i=is,ie
 
  174         sgn = 1 ; 
if (array(i,j)<0.0) sgn = -1
 
  177           ival = int(rs*i_pr(n), 8)
 
  179           ints_sum(n) = ints_sum(n) + sgn*ival
 
  185     if (
present(err)) 
then 
  187       if (overflow_error) &
 
  191       if (err > 0) 
then ; 
do n=1,
numint ; ints_sum(n) = 0 ;
 enddo ;
 endif 
  194         call mpp_error(fatal, 
"NaN in input field of mpp_reproducing_sum(_2d), this indicates numerical instability")
 
  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))
 
  200       if (overflow_error) 
then 
  201         call mpp_error(fatal, 
"Overflow in mpp_reproducing_sum(_2d).")
 
  211     do j=js,je ; 
do i=is,ie
 
  212       rsum(1) = rsum(1) + array(i,j);
 
  217     if (
present(err)) 
then ; err = 0 ;
 endif 
  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 
  226           write(mesg, 
'(ES13.5)') sum
 
  227           call mpp_error(fatal,
"Repro_sum_2d: Overflow in real_to_ints conversion of "//trim(mesg))
 
  233   if (
present(efp_sum)) efp_sum%v(:) = ints_sum(:)
 
  236     write(mesg,
'("2d RS: ", ES24.16, 6 Z17.16)') sum, ints_sum(1:
numint)
 
  243                             overflow_check, err) 
result(sum)
 
  244   real(r4_kind), 
dimension(:,:), 
intent(in) :: array
 
  245   integer,        
optional,          
intent(in) :: isr, ier, jsr, jer
 
  247   logical,        
optional,          
intent(in) :: reproducing
 
  248   logical,        
optional,          
intent(in) :: overflow_check
 
  249   integer,        
optional,         
intent(out) :: err
 
  252   real(r8_kind) :: array_r8(size(array,1), size(array,2))
 
  257                             overflow_check, err), r4_kind)
 
  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
 
  274   integer,            
optional,              
intent(out) :: err
 
  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
 
  285     "mpp_reproducing_sum: Too many processors are being used for the value of "//&
 
  286     "prec.  Reduce prec to (2^63-1)/mpp_npes.")
 
  288   prec_error = (2_8**62 + (2_8**62 - 1)) / 
mpp_npes()
 
  291   is = 1 ; ie = 
size(array,1) ; js = 1 ; je = 
size(array,2) ; ke = 
size(array,3)
 
  292   if (
present(isr)) 
then 
  294       "Value of isr too small in mpp_reproducing_sum(_3d).")
 
  297   if (
present(ier)) 
then 
  299       "Value of ier too large in mpp_reproducing_sum(_3d).")
 
  302   if (
present(jsr)) 
then 
  304       "Value of jsr too small in mpp_reproducing_sum(_3d).")
 
  307   if (
present(jer)) 
then 
  309       "Value of jer too large in mpp_reproducing_sum(_3d).")
 
  312   jsz = je+1-js; isz = ie+1-is
 
  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).")
 
  318     overflow_error = .false. ; nan_error = .false. ; max_mag_term = 0.0
 
  321         do j=js,je ; 
do i=is,ie
 
  327       do k=1,ke ; 
do j=js,je
 
  334       do k=1,ke ; 
do j=js,je ; 
do i=is,ie
 
  337       enddo ;
 enddo ;
 enddo 
  339     if (
present(err)) 
then 
  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 
  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))
 
  352       if (overflow_error) 
call mpp_error(fatal, 
"Overflow in mpp_reproducing_sum(_3d).")
 
  364     if (
present(efp_sum)) 
then 
  366       do k=1,ke ; 
call increment_ints(efp_sum%v(:), ints_sums(:,k)) ;
 enddo 
  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)
 
  377     overflow_error = .false. ; nan_error = .false. ; max_mag_term = 0.0
 
  380         do j=js,je ; 
do i=is,ie
 
  386       do k=1,ke ; 
do j=js,je
 
  393       do k=1,ke ; 
do j=js,je ; 
do i=is,ie
 
  396       enddo ;
 enddo ;
 enddo 
  398     if (
present(err)) 
then 
  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 
  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))
 
  411       if (overflow_error) 
call mpp_error(fatal, 
"Overflow in mpp_reproducing_sum(_3d).")
 
  419     if (
present(efp_sum)) efp_sum%v(:) = ints_sum(:)
 
  422       write(mesg,
'("3d RS: ", ES24.16, 6 Z17.16)') sum, ints_sum(1:
numint)
 
  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
 
  438   character(len=80) :: mesg
 
  439   integer(i8_kind) :: ival, prec_err
 
  442   prec_err = 
prec ; 
if (
present(prec_error)) prec_err = prec_error
 
  444   if ((r >= 1e30) .eqv. (r < 1e30)) 
then ; nan_error = .true. ; 
return ;
 endif 
  446   sgn = 1 ; 
if (r<0.0) sgn = -1
 
  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))
 
  458     ival = int(rs*i_pr(i), 8)
 
  467   integer(i8_kind), 
dimension(NUMINT), 
intent(in) :: ints
 
  473   do i=1,
numint ; r = r + pr(i)*ints(i) ;
 enddo 
  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
 
  486     int_sum(i) = int_sum(i) + int2(i)
 
  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
 
  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.
 
  500     if (abs(int_sum(1)) > 
prec) overflow_error = .true.
 
  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
 
  514   integer(i8_kind) :: ival
 
  517   if ((r >= 1e30) .eqv. (r < 1e30)) 
then ; nan_error = .true. ; 
return ;
 endif 
  518   sgn = 1 ; 
if (r<0.0) sgn = -1
 
  520   if (rs > abs(max_mag_term)) max_mag_term = r
 
  523     ival = int(rs*i_pr(i), 8)
 
  525     int_sum(i) = int_sum(i) + sgn*ival
 
  532   integer(i8_kind), 
dimension(NUMINT), 
intent(inout) :: int_sum
 
  533   integer(i8_kind),                
intent(in)    :: prec_error
 
  535   integer :: i, num_carry
 
  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
 
  542   if (abs(int_sum(1)) > prec_error) 
then 
  543     overflow_error = .true.
 
  551   integer(i8_kind), 
dimension(NUMINT), 
intent(inout) :: int_sum
 
  554   integer :: i, num_carry
 
  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
 
  565     if (abs(int_sum(i)) > 0) 
then 
  566       if (int_sum(i) < 0) positive = .false.
 
  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
 
  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
 
  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
 
  590 subroutine mpp_reset_efp_overflow_error()
 
  591   overflow_error = .false.
 
  592 end subroutine mpp_reset_efp_overflow_error
 
  594 function mpp_efp_plus(EFP1, EFP2)
 
  601 end function mpp_efp_plus
 
  603 function mpp_efp_minus(EFP1, EFP2)
 
  608   do i=1,
numint ; mpp_efp_minus%v(i) = -1*efp2%v(i) ;
 enddo 
  611 end function mpp_efp_minus
 
  616 subroutine mpp_efp_assign(EFP1, EFP2)
 
  621   do i=1,
numint ; efp1%v(i) = efp2%v(i) ;
 enddo 
  622 end subroutine mpp_efp_assign
 
  624 function mpp_efp_to_real(EFP1)
 
  626   real(r8_kind) :: mpp_efp_to_real
 
  630 end function mpp_efp_to_real
 
  632 function mpp_efp_real_diff(EFP1, EFP2)
 
  634   real(r8_kind) :: mpp_efp_real_diff
 
  638   efp_diff = efp1 - efp2
 
  639   mpp_efp_real_diff = mpp_efp_to_real(efp_diff)
 
  641 end function mpp_efp_real_diff
 
  643 function mpp_real_to_efp(val, overflow)
 
  644   real(r8_kind),    
intent(in) :: val
 
  645   logical, 
optional, 
intent(inout) :: overflow
 
  649   character(len=80) :: mesg
 
  651   if (
present(overflow)) 
then 
  652     mpp_real_to_efp%v(:) = 
real_to_ints(val, overflow=overflow)
 
  657       write(mesg, 
'(ES13.5)') val
 
  658       call mpp_error(fatal,
"Overflow in mpp_real_to_efp conversion of "//trim(mesg))
 
  662 end function mpp_real_to_efp
 
  668   integer, 
intent(in) :: nval
 
  669   logical, 
dimension(:), 
optional, 
intent(out) :: errors
 
  671   integer(i8_kind), 
dimension(NUMINT,nval) :: ints
 
  672   integer(i8_kind) :: prec_error
 
  673   logical :: error_found
 
  674   character(len=256) :: mesg
 
  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.")
 
  681   prec_error = (2_8**62 + (2_8**62 - 1)) / 
mpp_npes()
 
  683   overflow_error = .false. ; error_found = .false.
 
  685   do i=1,nval ; 
do n=1,
numint ; ints(n,i) = efps(i)%v(n) ;
 enddo ;
 enddo 
  689   if (
present(errors)) errors(:) = .false.
 
  691     overflow_error = .false.
 
  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)
 
  700     error_found = error_found .or. overflow_error
 
  702   if (error_found .and. .not.(
present(errors))) 
then 
  703     call mpp_error(fatal, 
"Overflow in mpp_efp_list_sum_across_PEs.")
 
  708 end module mpp_efp_mod
 
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...
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...
logical debug
Making this true enables debugging output.
subroutine carry_overflow(int_sum, prec_error)
This subroutine handles carrying of the overflow.
real(r8_kind), parameter i_prec
The inverse of prec.
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.
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.
real(r8_kind), parameter r_prec
A real version of prec.
integer, parameter numint
The number of long integers to use to represent a real number.
integer(i8_kind), parameter prec
The precision of each integer.
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.
subroutine increment_ints(int_sum, int2, prec_error)
This subroutine increments a number with another, both using the integer representation in real_to_in...
real(r4_kind) function mpp_reproducing_sum_r4_2d(array, isr, ier, jsr, jer, EFP_sum, reproducing, overflow_check, err)
real(r8_kind) function ints_to_real(ints)
This function reverses the conversion in real_to_ints.
subroutine regularize_ints(int_sum)
This subroutine carries the overflow, and then makes sure that all integers are of the same sign as t...
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,...
This interface uses a conversion to an integer representation of real numbers to give order-invariant...
The Extended Fixed Point (mpp_efp) type provides a public interface for doing sums and taking differe...
integer function mpp_npes()
Returns processor count for current pelist.
integer function mpp_pe()
Returns processor ID.