8 use fms_mod,
only : fms_end, mom_infra_init => fms_init
9 use memutils_mod,
only : print_memuse_stats
10 use mpp_mod,
only : pe_here => mpp_pe, root_pe => mpp_root_pe, num_pes => mpp_npes
11 use mpp_mod,
only : set_pelist => mpp_set_current_pelist, get_pelist => mpp_get_current_pelist
12 use mpp_mod,
only : broadcast => mpp_broadcast
13 use mpp_mod,
only : sum_across_pes => mpp_sum, min_across_pes => mpp_min
14 use mpp_mod,
only : max_across_pes => mpp_max
16 implicit none ;
private
18 public :: pe_here, root_pe, num_pes, mom_infra_init, mom_infra_end
19 public :: broadcast, sum_across_pes, min_across_pes, max_across_pes
21 public :: efp_plus, efp_minus, efp_to_real, real_to_efp, efp_real_diff
22 public ::
operator(+),
operator(-),
assignment(=)
23 public :: query_efp_overflow_error, reset_efp_overflow_error
24 public :: set_pelist, get_pelist
28 integer(kind=8),
parameter :: prec=2_8**46
29 real,
parameter :: r_prec=2.0**46
30 real,
parameter :: i_prec=1.0/(2.0**46)
31 integer,
parameter :: max_count_prec=2**(63-46)-1
36 integer,
parameter :: ni=6
38 real,
parameter,
dimension(ni) :: &
39 pr = (/ r_prec**2, r_prec, 1.0, 1.0/r_prec, 1.0/r_prec**2, 1.0/r_prec**3 /)
41 real,
parameter,
dimension(ni) :: &
42 i_pr = (/ 1.0/r_prec**2, 1.0/r_prec, 1.0, r_prec, r_prec**2, r_prec**3 /)
44 real,
parameter :: max_efp_float = pr(1) * (2.**63 - 1.)
49 logical :: overflow_error = .false.
50 logical :: nan_error = .false.
51 logical :: debug = .false.
55 module procedure reproducing_sum_2d, reproducing_sum_3d
65 integer(kind=8),
dimension(ni) :: v
69 interface operator (+) ; module
procedure EFP_plus ; end interface
71 interface operator (-) ; module
procedure EFP_minus ; end interface
73 interface assignment(=); module
procedure EFP_assign ; end interface
81 function reproducing_sum_2d(array, isr, ier, jsr, jer, EFP_sum, reproducing, &
82 overflow_check, err)
result(sum)
83 real,
dimension(:,:),
intent(in) :: array
84 integer,
optional,
intent(in) :: isr
86 integer,
optional,
intent(in) :: ier
88 integer,
optional,
intent(in) :: jsr
90 integer,
optional,
intent(in) :: jer
92 type(
efp_type),
optional,
intent(out) :: efp_sum
93 logical,
optional,
intent(in) :: reproducing
95 logical,
optional,
intent(in) :: overflow_check
99 integer,
optional,
intent(out) :: err
108 integer(kind=8),
dimension(ni) :: ints_sum
109 integer(kind=8) :: ival, prec_error
112 logical :: repro, over_check
113 character(len=256) :: mesg
114 integer :: i, j, n, is, ie, js, je, sgn
116 if (num_pes() > max_count_prec)
call mom_error(fatal, &
117 "reproducing_sum: Too many processors are being used for the value of "//&
118 "prec. Reduce prec to (2^63-1)/num_PEs.")
120 prec_error = (2_8**62 + (2_8**62 - 1)) / num_pes()
122 is = 1 ; ie =
size(array,1) ; js = 1 ; je =
size(array,2 )
123 if (
present(isr))
then
124 if (isr < is)
call mom_error(fatal, &
125 "Value of isr too small in reproducing_sum_2d.")
128 if (
present(ier))
then
129 if (ier > ie)
call mom_error(fatal, &
130 "Value of ier too large in reproducing_sum_2d.")
133 if (
present(jsr))
then
134 if (jsr < js)
call mom_error(fatal, &
135 "Value of jsr too small in reproducing_sum_2d.")
138 if (
present(jer))
then
139 if (jer > je)
call mom_error(fatal, &
140 "Value of jer too large in reproducing_sum_2d.")
144 repro = .true. ;
if (
present(reproducing)) repro = reproducing
145 over_check = .true. ;
if (
present(overflow_check)) over_check = overflow_check
148 overflow_error = .false. ; nan_error = .false. ; max_mag_term = 0.0
151 if ((je+1-js)*(ie+1-is) < max_count_prec)
then
152 do j=js,je ;
do i=is,ie
153 call increment_ints_faster(ints_sum, array(i,j), max_mag_term)
155 call carry_overflow(ints_sum, prec_error)
156 elseif ((ie+1-is) < max_count_prec)
then
159 call increment_ints_faster(ints_sum, array(i,j), max_mag_term)
161 call carry_overflow(ints_sum, prec_error)
164 do j=js,je ;
do i=is,ie
165 call increment_ints(ints_sum, real_to_ints(array(i,j), prec_error), &
170 do j=js,je ;
do i=is,ie
171 sgn = 1 ;
if (array(i,j)<0.0) sgn = -1
174 ival = int(rs*i_pr(n), 8)
176 ints_sum(n) = ints_sum(n) + sgn*ival
179 call carry_overflow(ints_sum, prec_error)
182 if (
present(err))
then
184 if (overflow_error) &
188 if (err > 0)
then ;
do n=1,ni ; ints_sum(n) = 0 ;
enddo ;
endif
191 call mom_error(fatal,
"NaN in input field of reproducing_sum(_2d).")
193 if (abs(max_mag_term) >= prec_error*pr(1))
then
194 write(mesg,
'(ES13.5)') max_mag_term
195 call mom_error(fatal,
"Overflow in reproducing_sum(_2d) conversion of "//trim(mesg))
197 if (overflow_error)
then
198 call mom_error(fatal,
"Overflow in reproducing_sum(_2d).")
202 call sum_across_pes(ints_sum, ni)
204 call regularize_ints(ints_sum)
205 sum = ints_to_real(ints_sum)
208 do j=js,je ;
do i=is,ie
209 rsum(1) = rsum(1) + array(i,j)
211 call sum_across_pes(rsum,1)
214 if (
present(err))
then ; err = 0 ;
endif
216 if (debug .or.
present(efp_sum))
then
217 overflow_error = .false.
218 ints_sum = real_to_ints(sum, prec_error, overflow_error)
219 if (overflow_error)
then
220 if (
present(err))
then
223 write(mesg,
'(ES13.5)') sum
224 call mom_error(fatal,
"Repro_sum_2d: Overflow in real_to_ints conversion of "//trim(mesg))
230 if (
present(efp_sum)) efp_sum%v(:) = ints_sum(:)
233 write(mesg,
'("2d RS: ", ES24.16, 6 Z17.16)') sum, ints_sum(1:ni)
234 call mom_mesg(mesg, 3)
237 end function reproducing_sum_2d
243 function reproducing_sum_3d(array, isr, ier, jsr, jer, sums, EFP_sum, err) &
245 real,
dimension(:,:,:),
intent(in) :: array
246 integer,
optional,
intent(in) :: isr
248 integer,
optional,
intent(in) :: ier
250 integer,
optional,
intent(in) :: jsr
252 integer,
optional,
intent(in) :: jer
254 real,
dimension(:),
optional,
intent(out) :: sums
255 type(
efp_type),
optional,
intent(out) :: efp_sum
256 integer,
optional,
intent(out) :: err
266 integer(kind=8),
dimension(ni) :: ints_sum
267 integer(kind=8),
dimension(ni,size(array,3)) :: ints_sums
268 integer(kind=8) :: prec_error
269 character(len=256) :: mesg
270 integer :: i, j, k, is, ie, js, je, ke, isz, jsz, n
272 if (num_pes() > max_count_prec)
call mom_error(fatal, &
273 "reproducing_sum: Too many processors are being used for the value of "//&
274 "prec. Reduce prec to (2^63-1)/num_PEs.")
276 prec_error = (2_8**62 + (2_8**62 - 1)) / num_pes()
279 is = 1 ; ie =
size(array,1) ; js = 1 ; je =
size(array,2) ; ke =
size(array,3)
280 if (
present(isr))
then
281 if (isr < is)
call mom_error(fatal, &
282 "Value of isr too small in reproducing_sum(_3d).")
285 if (
present(ier))
then
286 if (ier > ie)
call mom_error(fatal, &
287 "Value of ier too large in reproducing_sum(_3d).")
290 if (
present(jsr))
then
291 if (jsr < js)
call mom_error(fatal, &
292 "Value of jsr too small in reproducing_sum(_3d).")
295 if (
present(jer))
then
296 if (jer > je)
call mom_error(fatal, &
297 "Value of jer too large in reproducing_sum(_3d).")
300 jsz = je+1-js; isz = ie+1-is
302 if (
present(sums))
then
303 if (
size(sums) > ke)
call mom_error(fatal,
"Sums is smaller than "//&
304 "the vertical extent of array in reproducing_sum(_3d).")
306 overflow_error = .false. ; nan_error = .false. ; max_mag_term = 0.0
307 if (jsz*isz < max_count_prec)
then
309 do j=js,je ;
do i=is,ie
310 call increment_ints_faster(ints_sums(:,k), array(i,j,k), max_mag_term)
312 call carry_overflow(ints_sums(:,k), prec_error)
314 elseif (isz < max_count_prec)
then
315 do k=1,ke ;
do j=js,je
317 call increment_ints_faster(ints_sums(:,k), array(i,j,k), max_mag_term)
319 call carry_overflow(ints_sums(:,k), prec_error)
322 do k=1,ke ;
do j=js,je ;
do i=is,ie
323 call increment_ints(ints_sums(:,k), &
324 real_to_ints(array(i,j,k), prec_error), prec_error)
325 enddo ;
enddo ;
enddo
327 if (
present(err))
then
329 if (abs(max_mag_term) >= prec_error*pr(1)) err = err+1
330 if (overflow_error) err = err+2
331 if (nan_error) err = err+2
332 if (err > 0)
then ;
do k=1,ke ;
do n=1,ni ; ints_sums(n,k) = 0 ;
enddo ;
enddo ;
endif
334 if (nan_error)
call mom_error(fatal,
"NaN in input field of reproducing_sum(_3d).")
335 if (abs(max_mag_term) >= prec_error*pr(1))
then
336 write(mesg,
'(ES13.5)') max_mag_term
337 call mom_error(fatal,
"Overflow in reproducing_sum(_3d) conversion of "//trim(mesg))
339 if (overflow_error)
call mom_error(fatal,
"Overflow in reproducing_sum(_3d).")
342 call sum_across_pes(ints_sums(:,1:ke), ni*ke)
346 call regularize_ints(ints_sums(:,k))
347 sums(k) = ints_to_real(ints_sums(:,k))
351 if (
present(efp_sum))
then
353 do k=1,ke ;
call increment_ints(efp_sum%v(:), ints_sums(:,k)) ;
enddo
357 do n=1,ni ; ints_sum(n) = 0 ;
enddo
358 do k=1,ke ;
do n=1,ni ; ints_sum(n) = ints_sum(n) + ints_sums(n,k) ;
enddo ;
enddo
359 write(mesg,
'("3D RS: ", ES24.16, 6 Z17.16)') sum, ints_sum(1:ni)
360 call mom_mesg(mesg, 3)
364 overflow_error = .false. ; nan_error = .false. ; max_mag_term = 0.0
365 if (jsz*isz < max_count_prec)
then
367 do j=js,je ;
do i=is,ie
368 call increment_ints_faster(ints_sum, array(i,j,k), max_mag_term)
370 call carry_overflow(ints_sum, prec_error)
372 elseif (isz < max_count_prec)
then
373 do k=1,ke ;
do j=js,je
375 call increment_ints_faster(ints_sum, array(i,j,k), max_mag_term)
377 call carry_overflow(ints_sum, prec_error)
380 do k=1,ke ;
do j=js,je ;
do i=is,ie
381 call increment_ints(ints_sum, real_to_ints(array(i,j,k), prec_error), &
383 enddo ;
enddo ;
enddo
385 if (
present(err))
then
387 if (abs(max_mag_term) >= prec_error*pr(1)) err = err+1
388 if (overflow_error) err = err+2
389 if (nan_error) err = err+2
390 if (err > 0)
then ;
do n=1,ni ; ints_sum(n) = 0 ;
enddo ;
endif
392 if (nan_error)
call mom_error(fatal,
"NaN in input field of reproducing_sum(_3d).")
393 if (abs(max_mag_term) >= prec_error*pr(1))
then
394 write(mesg,
'(ES13.5)') max_mag_term
395 call mom_error(fatal,
"Overflow in reproducing_sum(_3d) conversion of "//trim(mesg))
397 if (overflow_error)
call mom_error(fatal,
"Overflow in reproducing_sum(_3d).")
400 call sum_across_pes(ints_sum, ni)
402 call regularize_ints(ints_sum)
403 sum = ints_to_real(ints_sum)
405 if (
present(efp_sum)) efp_sum%v(:) = ints_sum(:)
408 write(mesg,
'("3d RS: ", ES24.16, 6 Z17.16)') sum, ints_sum(1:ni)
409 call mom_mesg(mesg, 3)
413 end function reproducing_sum_3d
416 function real_to_ints(r, prec_error, overflow)
result(ints)
417 real,
intent(in) :: r
418 integer(kind=8),
optional,
intent(in) :: prec_error
422 logical,
optional,
intent(inout) :: overflow
424 integer(kind=8),
dimension(ni) :: ints
429 character(len=80) :: mesg
430 integer(kind=8) :: ival, prec_err
433 prec_err = prec ;
if (
present(prec_error)) prec_err = prec_error
435 if ((r >= 1e30) .eqv. (r < 1e30))
then ; nan_error = .true. ;
return ;
endif
437 sgn = 1 ;
if (r<0.0) sgn = -1
440 if (
present(overflow))
then
441 if (.not.(rs < prec_err*pr(1))) overflow = .true.
442 if ((r >= 1e30) .eqv. (r < 1e30)) overflow = .true.
443 elseif (.not.(rs < prec_err*pr(1)))
then
444 write(mesg,
'(ES13.5)') r
445 call mom_error(fatal,
"Overflow in real_to_ints conversion of "//trim(mesg))
449 ival = int(rs*i_pr(i), 8)
454 end function real_to_ints
458 function ints_to_real(ints)
result(r)
459 integer(kind=8),
dimension(ni),
intent(in) :: ints
466 do i=1,ni ; r = r + pr(i)*ints(i) ;
enddo
467 end function ints_to_real
471 subroutine increment_ints(int_sum, int2, prec_error)
472 integer(kind=8),
dimension(ni),
intent(inout) :: int_sum
473 integer(kind=8),
dimension(ni),
intent(in) :: int2
474 integer(kind=8),
optional,
intent(in) :: prec_error
484 int_sum(i) = int_sum(i) + int2(i)
486 if (int_sum(i) > prec)
then
487 int_sum(i) = int_sum(i) - prec
488 int_sum(i-1) = int_sum(i-1) + 1
489 elseif (int_sum(i) < -prec)
then
490 int_sum(i) = int_sum(i) + prec
491 int_sum(i-1) = int_sum(i-1) - 1
494 int_sum(1) = int_sum(1) + int2(1)
495 if (
present(prec_error))
then
496 if (abs(int_sum(1)) > prec_error) overflow_error = .true.
498 if (abs(int_sum(1)) > prec) overflow_error = .true.
501 end subroutine increment_ints
505 subroutine increment_ints_faster(int_sum, r, max_mag_term)
506 integer(kind=8),
dimension(ni),
intent(inout) :: int_sum
507 real,
intent(in) :: r
508 real,
intent(inout) :: max_mag_term
514 integer(kind=8) :: 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 if (rs > max_efp_float)
then
524 overflow_error = .true.
529 ival = int(rs*i_pr(i), 8)
531 int_sum(i) = int_sum(i) + sgn*ival
534 end subroutine increment_ints_faster
537 subroutine carry_overflow(int_sum, prec_error)
538 integer(kind=8),
dimension(ni),
intent(inout) :: int_sum
540 integer(kind=8),
intent(in) :: prec_error
546 integer :: i, num_carry
548 do i=ni,2,-1 ;
if (abs(int_sum(i)) >= prec)
then
549 num_carry = int(int_sum(i) * i_prec)
550 int_sum(i) = int_sum(i) - num_carry*prec
551 int_sum(i-1) = int_sum(i-1) + num_carry
553 if (abs(int_sum(1)) > prec_error)
then
554 overflow_error = .true.
557 end subroutine carry_overflow
561 subroutine regularize_ints(int_sum)
562 integer(kind=8),
dimension(ni), &
563 intent(inout) :: int_sum
570 integer :: i, num_carry
572 do i=ni,2,-1 ;
if (abs(int_sum(i)) >= prec)
then
573 num_carry = int(int_sum(i) * i_prec)
574 int_sum(i) = int_sum(i) - num_carry*prec
575 int_sum(i-1) = int_sum(i-1) + num_carry
581 if (abs(int_sum(i)) > 0)
then
582 if (int_sum(i) < 0) positive = .false.
588 do i=ni,2,-1 ;
if (int_sum(i) < 0)
then
589 int_sum(i) = int_sum(i) + prec
590 int_sum(i-1) = int_sum(i-1) - 1
593 do i=ni,2,-1 ;
if (int_sum(i) > 0)
then
594 int_sum(i) = int_sum(i) - prec
595 int_sum(i-1) = int_sum(i-1) + 1
599 end subroutine regularize_ints
602 function query_efp_overflow_error()
603 logical :: query_efp_overflow_error
604 query_efp_overflow_error = overflow_error
605 end function query_efp_overflow_error
608 subroutine reset_efp_overflow_error()
609 overflow_error = .false.
610 end subroutine reset_efp_overflow_error
613 function efp_plus(EFP1, EFP2)
620 call increment_ints(efp_plus%v(:), efp2%v(:))
621 end function efp_plus
624 function efp_minus(EFP1, EFP2)
631 do i=1,ni ; efp_minus%v(i) = -1*efp2%v(i) ;
enddo
633 call increment_ints(efp_minus%v(:), efp1%v(:))
634 end function efp_minus
637 subroutine efp_assign(EFP1, EFP2)
645 do i=1,ni ; efp1%v(i) = efp2%v(i) ;
enddo
646 end subroutine efp_assign
649 function efp_to_real(EFP1)
653 call regularize_ints(efp1%v)
654 efp_to_real = ints_to_real(efp1%v)
655 end function efp_to_real
659 function efp_real_diff(EFP1, EFP2)
663 real :: efp_real_diff
667 efp_diff = efp1 - efp2
668 efp_real_diff = efp_to_real(efp_diff)
670 end function efp_real_diff
673 function real_to_efp(val, overflow)
674 real,
intent(in) :: val
675 logical,
optional,
intent(inout) :: overflow
680 character(len=80) :: mesg
682 if (
present(overflow))
then
683 real_to_efp%v(:) = real_to_ints(val, overflow=overflow)
686 real_to_efp%v(:) = real_to_ints(val, overflow=over)
688 write(mesg,
'(ES13.5)') val
689 call mom_error(fatal,
"Overflow in real_to_EFP conversion of "//trim(mesg))
693 end function real_to_efp
697 subroutine efp_list_sum_across_pes(EFPs, nval, errors)
699 intent(inout) :: efps
701 integer,
intent(in) :: nval
702 logical,
dimension(:), &
703 optional,
intent(out) :: errors
708 integer(kind=8),
dimension(ni,nval) :: ints
709 integer(kind=8) :: prec_error
710 logical :: error_found
711 character(len=256) :: mesg
714 if (num_pes() > max_count_prec)
call mom_error(fatal, &
715 "reproducing_sum: Too many processors are being used for the value of "//&
716 "prec. Reduce prec to (2^63-1)/num_PEs.")
718 prec_error = (2_8**62 + (2_8**62 - 1)) / num_pes()
720 overflow_error = .false. ; error_found = .false.
722 do i=1,nval ;
do n=1,ni ; ints(n,i) = efps(i)%v(n) ;
enddo ;
enddo
724 call sum_across_pes(ints(:,:), ni*nval)
726 if (
present(errors)) errors(:) = .false.
728 overflow_error = .false.
729 call carry_overflow(ints(:,i), prec_error)
730 do n=1,ni ; efps(i)%v(n) = ints(n,i) ;
enddo
731 if (
present(errors)) errors(i) = overflow_error
732 if (overflow_error)
then
733 write (mesg,
'("EFP_list_sum_across_PEs error at ",i6," val was ",ES12.6, ", prec_error = ",ES12.6)') &
734 i, efp_to_real(efps(i)), real(prec_error)
735 call mom_error(warning, mesg)
737 error_found = error_found .or. overflow_error
739 if (error_found .and. .not.(
present(errors)))
then
740 call mom_error(fatal,
"Overflow in EFP_list_sum_across_PEs.")
743 end subroutine efp_list_sum_across_pes
747 subroutine mom_infra_end
748 call print_memuse_stats(
'Memory HiWaterMark', always=.true. )
750 end subroutine mom_infra_end