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
22 public ::
operator(+),
operator(-),
assignment(=)
24 public :: set_pelist, get_pelist
28 integer(kind=8),
parameter ::
prec=2_8**46
30 real,
parameter ::
i_prec=1.0/(2.0**46)
36 integer,
parameter ::
ni=6
38 real,
parameter,
dimension(ni) :: &
41 real,
parameter,
dimension(ni) :: &
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
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
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
125 "Value of isr too small in reproducing_sum_2d.")
128 if (
present(ier))
then
130 "Value of ier too large in reproducing_sum_2d.")
133 if (
present(jsr))
then
135 "Value of jsr too small in reproducing_sum_2d.")
138 if (
present(jer))
then
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
152 do j=js,je ;
do i=is,ie
164 do j=js,je ;
do i=is,ie
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
182 if (
present(err))
then
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))
198 call mom_error(fatal,
"Overflow in reproducing_sum(_2d).")
202 call sum_across_pes(ints_sum,
ni)
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
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)
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
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
282 "Value of isr too small in reproducing_sum(_3d).")
285 if (
present(ier))
then
287 "Value of ier too large in reproducing_sum(_3d).")
290 if (
present(jsr))
then
292 "Value of jsr too small in reproducing_sum(_3d).")
295 if (
present(jer))
then
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).")
309 do j=js,je ;
do i=is,ie
315 do k=1,ke ;
do j=js,je
322 do k=1,ke ;
do j=js,je ;
do i=is,ie
325 enddo ;
enddo ;
enddo
327 if (
present(err))
then
329 if (abs(max_mag_term) >= prec_error*
pr(1)) err = err+1
332 if (err > 0)
then ;
do k=1,ke ;
do n=1,
ni ; ints_sums(n,k) = 0 ;
enddo ;
enddo ;
endif
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))
342 call sum_across_pes(ints_sums(:,1:ke),
ni*ke)
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)
367 do j=js,je ;
do i=is,ie
373 do k=1,ke ;
do j=js,je
380 do k=1,ke ;
do j=js,je ;
do i=is,ie
383 enddo ;
enddo ;
enddo
385 if (
present(err))
then
387 if (abs(max_mag_term) >= prec_error*
pr(1)) err = err+1
390 if (err > 0)
then ;
do n=1,
ni ; ints_sum(n) = 0 ;
enddo ;
endif
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))
400 call sum_across_pes(ints_sum,
ni)
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)
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)
459 integer(kind=8),
dimension(ni),
intent(in) :: ints
466 do i=1,
ni ; r = r +
pr(i)*ints(i) ;
enddo
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
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
529 ival = int(rs*
i_pr(i), 8)
531 int_sum(i) = int_sum(i) + sgn*ival
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
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
613 function efp_plus(EFP1, EFP2)
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
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
667 efp_diff = efp1 - efp2
674 real,
intent(in) :: val
675 logical,
optional,
intent(inout) :: overflow
680 character(len=80) :: mesg
682 if (
present(overflow))
then
688 write(mesg,
'(ES13.5)') val
689 call mom_error(fatal,
"Overflow in real_to_EFP conversion of "//trim(mesg))
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
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()
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.
730 do n=1,
ni ; efps(i)%v(n) = ints(n,i) ;
enddo
733 write (mesg,
'("EFP_list_sum_across_PEs error at ",i6," val was ",ES12.6, ", prec_error = ",ES12.6)') &
739 if (error_found .and. .not.(
present(errors)))
then
740 call mom_error(fatal,
"Overflow in EFP_list_sum_across_PEs.")
748 call print_memuse_stats(
'Memory HiWaterMark', always=.true. )