MOM6
MOM_coms.F90
Go to the documentation of this file.
1 !> Interfaces to non-domain-oriented communication subroutines, including the
2 !! MOM6 reproducing sums facility
3 module mom_coms
4 
5 ! This file is part of MOM6. See LICENSE.md for the license.
6 
7 use mom_error_handler, only : mom_error, mom_mesg, fatal, warning
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
15 
16 implicit none ; private
17 
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
25 ! This module provides interfaces to the non-domain-oriented communication
26 ! subroutines.
27 
28 integer(kind=8), parameter :: prec=2_8**46 !< The precision of each integer.
29 real, parameter :: r_prec=2.0**46 !< A real version of prec.
30 real, parameter :: i_prec=1.0/(2.0**46) !< The inverse of prec.
31 integer, parameter :: max_count_prec=2**(63-46)-1
32  !< The number of values that can be added together
33  !! with the current value of prec before there will
34  !! be roundoff problems.
35 
36 integer, parameter :: ni=6 !< The number of long integers to use to represent
37  !< a real number.
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 /)
40  !< An array of the real precision of each of the integers
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 /)
43  !< An array of the inverse of the real precision of each of the integers
44 real, parameter :: max_efp_float = pr(1) * (2.**63 - 1.)
45  !< The largest float with an EFP representation.
46  !! NOTE: Only the first bin can exceed precision,
47  !! but is bounded by the largest signed integer.
48 
49 logical :: overflow_error = .false. !< This becomes true if an overflow is encountered.
50 logical :: nan_error = .false. !< This becomes true if a NaN is encountered.
51 logical :: debug = .false. !< Making this true enables debugging output.
52 
53 !> Find an accurate and order-invariant sum of distributed 2d or 3d fields
54 interface reproducing_sum
55  module procedure reproducing_sum_2d, reproducing_sum_3d
56 end interface reproducing_sum
57 
58 !> The Extended Fixed Point (EFP) type provides a public interface for doing sums
59 !! and taking differences with this type.
60 !!
61 !! The use of this type is documented in
62 !! Hallberg, R. & A. Adcroft, 2014: An Order-invariant Real-to-Integer Conversion Sum.
63 !! Parallel Computing, 40(5-6), doi:10.1016/j.parco.2014.04.007.
64 type, public :: efp_type ; private
65  integer(kind=8), dimension(ni) :: v !< The value in this type
66 end type efp_type
67 
68 !> Add two extended-fixed-point numbers
69 interface operator (+) ; module procedure EFP_plus ; end interface
70 !> Subtract one extended-fixed-point number from another
71 interface operator (-) ; module procedure EFP_minus ; end interface
72 !> Copy the value of one extended-fixed-point number into another
73 interface assignment(=); module procedure EFP_assign ; end interface
74 
75 contains
76 
77 !> This subroutine uses a conversion to an integer representation of real numbers to give an
78 !! order-invariant sum of distributed 2-D arrays that reproduces across domain decomposition.
79 !! This technique is described in Hallberg & Adcroft, 2014, Parallel Computing,
80 !! doi:10.1016/j.parco.2014.04.007.
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 !< The array to be summed
84  integer, optional, intent(in) :: isr !< The starting i-index of the sum, noting
85  !! that the array indices starts at 1
86  integer, optional, intent(in) :: ier !< The ending i-index of the sum, noting
87  !! that the array indices starts at 1
88  integer, optional, intent(in) :: jsr !< The starting j-index of the sum, noting
89  !! that the array indices starts at 1
90  integer, optional, intent(in) :: jer !< The ending j-index of the sum, noting
91  !! that the array indices starts at 1
92  type(efp_type), optional, intent(out) :: efp_sum !< The result in extended fixed point format
93  logical, optional, intent(in) :: reproducing !< If present and false, do the sum
94  !! using the naive non-reproducing approach
95  logical, optional, intent(in) :: overflow_check !< If present and false, disable
96  !! checking for overflows in incremental results.
97  !! This can speed up calculations if the number
98  !! of values being summed is small enough
99  integer, optional, intent(out) :: err !< If present, return an error code instead of
100  !! triggering any fatal errors directly from
101  !! this routine.
102  real :: sum !< Result
103 
104  ! This subroutine uses a conversion to an integer representation
105  ! of real numbers to give order-invariant sums that will reproduce
106  ! across PE count. This idea comes from R. Hallberg and A. Adcroft.
107 
108  integer(kind=8), dimension(ni) :: ints_sum
109  integer(kind=8) :: ival, prec_error
110  real :: rsum(1), rs
111  real :: max_mag_term
112  logical :: repro, over_check
113  character(len=256) :: mesg
114  integer :: i, j, n, is, ie, js, je, sgn
115 
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.")
119 
120  prec_error = (2_8**62 + (2_8**62 - 1)) / num_pes()
121 
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.")
126  is = isr
127  endif
128  if (present(ier)) then
129  if (ier > ie) call mom_error(fatal, &
130  "Value of ier too large in reproducing_sum_2d.")
131  ie = ier
132  endif
133  if (present(jsr)) then
134  if (jsr < js) call mom_error(fatal, &
135  "Value of jsr too small in reproducing_sum_2d.")
136  js = jsr
137  endif
138  if (present(jer)) then
139  if (jer > je) call mom_error(fatal, &
140  "Value of jer too large in reproducing_sum_2d.")
141  je = jer
142  endif
143 
144  repro = .true. ; if (present(reproducing)) repro = reproducing
145  over_check = .true. ; if (present(overflow_check)) over_check = overflow_check
146 
147  if (repro) then
148  overflow_error = .false. ; nan_error = .false. ; max_mag_term = 0.0
149  ints_sum(:) = 0
150  if (over_check) then
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)
154  enddo ; enddo
155  call carry_overflow(ints_sum, prec_error)
156  elseif ((ie+1-is) < max_count_prec) then
157  do j=js,je
158  do i=is,ie
159  call increment_ints_faster(ints_sum, array(i,j), max_mag_term)
160  enddo
161  call carry_overflow(ints_sum, prec_error)
162  enddo
163  else
164  do j=js,je ; do i=is,ie
165  call increment_ints(ints_sum, real_to_ints(array(i,j), prec_error), &
166  prec_error)
167  enddo ; enddo
168  endif
169  else
170  do j=js,je ; do i=is,ie
171  sgn = 1 ; if (array(i,j)<0.0) sgn = -1
172  rs = abs(array(i,j))
173  do n=1,ni
174  ival = int(rs*i_pr(n), 8)
175  rs = rs - ival*pr(n)
176  ints_sum(n) = ints_sum(n) + sgn*ival
177  enddo
178  enddo ; enddo
179  call carry_overflow(ints_sum, prec_error)
180  endif
181 
182  if (present(err)) then
183  err = 0
184  if (overflow_error) &
185  err = err+2
186  if (nan_error) &
187  err = err+4
188  if (err > 0) then ; do n=1,ni ; ints_sum(n) = 0 ; enddo ; endif
189  else
190  if (nan_error) then
191  call mom_error(fatal, "NaN in input field of reproducing_sum(_2d).")
192  endif
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))
196  endif
197  if (overflow_error) then
198  call mom_error(fatal, "Overflow in reproducing_sum(_2d).")
199  endif
200  endif
201 
202  call sum_across_pes(ints_sum, ni)
203 
204  call regularize_ints(ints_sum)
205  sum = ints_to_real(ints_sum)
206  else
207  rsum(1) = 0.0
208  do j=js,je ; do i=is,ie
209  rsum(1) = rsum(1) + array(i,j)
210  enddo ; enddo
211  call sum_across_pes(rsum,1)
212  sum = rsum(1)
213 
214  if (present(err)) then ; err = 0 ; endif
215 
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
221  err = err + 2
222  else
223  write(mesg, '(ES13.5)') sum
224  call mom_error(fatal,"Repro_sum_2d: Overflow in real_to_ints conversion of "//trim(mesg))
225  endif
226  endif
227  endif
228  endif
229 
230  if (present(efp_sum)) efp_sum%v(:) = ints_sum(:)
231 
232  if (debug) then
233  write(mesg,'("2d RS: ", ES24.16, 6 Z17.16)') sum, ints_sum(1:ni)
234  call mom_mesg(mesg, 3)
235  endif
236 
237 end function reproducing_sum_2d
238 
239 !> This subroutine uses a conversion to an integer representation of real numbers to give an
240 !! order-invariant sum of distributed 3-D arrays that reproduces across domain decomposition.
241 !! This technique is described in Hallberg & Adcroft, 2014, Parallel Computing,
242 !! doi:10.1016/j.parco.2014.04.007.
243 function reproducing_sum_3d(array, isr, ier, jsr, jer, sums, EFP_sum, err) &
244  result(sum)
245  real, dimension(:,:,:), intent(in) :: array !< The array to be summed
246  integer, optional, intent(in) :: isr !< The starting i-index of the sum, noting
247  !! that the array indices starts at 1
248  integer, optional, intent(in) :: ier !< The ending i-index of the sum, noting
249  !! that the array indices starts at 1
250  integer, optional, intent(in) :: jsr !< The starting j-index of the sum, noting
251  !! that the array indices starts at 1
252  integer, optional, intent(in) :: jer !< The ending j-index of the sum, noting
253  !! that the array indices starts at 1
254  real, dimension(:), optional, intent(out) :: sums !< The sums by vertical layer
255  type(efp_type), optional, intent(out) :: efp_sum !< The result in extended fixed point format
256  integer, optional, intent(out) :: err !< If present, return an error code instead of
257  !! triggering any fatal errors directly from
258  !! this routine.
259  real :: sum !< Result
260 
261  ! This subroutine uses a conversion to an integer representation
262  ! of real numbers to give order-invariant sums that will reproduce
263  ! across PE count. This idea comes from R. Hallberg and A. Adcroft.
264 
265  real :: max_mag_term
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
271 
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.")
275 
276  prec_error = (2_8**62 + (2_8**62 - 1)) / num_pes()
277  max_mag_term = 0.0
278 
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).")
283  is = isr
284  endif
285  if (present(ier)) then
286  if (ier > ie) call mom_error(fatal, &
287  "Value of ier too large in reproducing_sum(_3d).")
288  ie = ier
289  endif
290  if (present(jsr)) then
291  if (jsr < js) call mom_error(fatal, &
292  "Value of jsr too small in reproducing_sum(_3d).")
293  js = jsr
294  endif
295  if (present(jer)) then
296  if (jer > je) call mom_error(fatal, &
297  "Value of jer too large in reproducing_sum(_3d).")
298  je = jer
299  endif
300  jsz = je+1-js; isz = ie+1-is
301 
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).")
305  ints_sums(:,:) = 0
306  overflow_error = .false. ; nan_error = .false. ; max_mag_term = 0.0
307  if (jsz*isz < max_count_prec) then
308  do k=1,ke
309  do j=js,je ; do i=is,ie
310  call increment_ints_faster(ints_sums(:,k), array(i,j,k), max_mag_term)
311  enddo ; enddo
312  call carry_overflow(ints_sums(:,k), prec_error)
313  enddo
314  elseif (isz < max_count_prec) then
315  do k=1,ke ; do j=js,je
316  do i=is,ie
317  call increment_ints_faster(ints_sums(:,k), array(i,j,k), max_mag_term)
318  enddo
319  call carry_overflow(ints_sums(:,k), prec_error)
320  enddo ; enddo
321  else
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
326  endif
327  if (present(err)) then
328  err = 0
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
333  else
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))
338  endif
339  if (overflow_error) call mom_error(fatal, "Overflow in reproducing_sum(_3d).")
340  endif
341 
342  call sum_across_pes(ints_sums(:,1:ke), ni*ke)
343 
344  sum = 0.0
345  do k=1,ke
346  call regularize_ints(ints_sums(:,k))
347  sums(k) = ints_to_real(ints_sums(:,k))
348  sum = sum + sums(k)
349  enddo
350 
351  if (present(efp_sum)) then
352  efp_sum%v(:) = 0
353  do k=1,ke ; call increment_ints(efp_sum%v(:), ints_sums(:,k)) ; enddo
354  endif
355 
356  if (debug) then
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)
361  endif
362  else
363  ints_sum(:) = 0
364  overflow_error = .false. ; nan_error = .false. ; max_mag_term = 0.0
365  if (jsz*isz < max_count_prec) then
366  do k=1,ke
367  do j=js,je ; do i=is,ie
368  call increment_ints_faster(ints_sum, array(i,j,k), max_mag_term)
369  enddo ; enddo
370  call carry_overflow(ints_sum, prec_error)
371  enddo
372  elseif (isz < max_count_prec) then
373  do k=1,ke ; do j=js,je
374  do i=is,ie
375  call increment_ints_faster(ints_sum, array(i,j,k), max_mag_term)
376  enddo
377  call carry_overflow(ints_sum, prec_error)
378  enddo ; enddo
379  else
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), &
382  prec_error)
383  enddo ; enddo ; enddo
384  endif
385  if (present(err)) then
386  err = 0
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
391  else
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))
396  endif
397  if (overflow_error) call mom_error(fatal, "Overflow in reproducing_sum(_3d).")
398  endif
399 
400  call sum_across_pes(ints_sum, ni)
401 
402  call regularize_ints(ints_sum)
403  sum = ints_to_real(ints_sum)
404 
405  if (present(efp_sum)) efp_sum%v(:) = ints_sum(:)
406 
407  if (debug) then
408  write(mesg,'("3d RS: ", ES24.16, 6 Z17.16)') sum, ints_sum(1:ni)
409  call mom_mesg(mesg, 3)
410  endif
411  endif
412 
413 end function reproducing_sum_3d
414 
415 !> Convert a real number into the array of integers constitute its extended-fixed-point representation
416 function real_to_ints(r, prec_error, overflow) result(ints)
417  real, intent(in) :: r !< The real number being converted
418  integer(kind=8), optional, intent(in) :: prec_error !< The PE-count dependent precision of the
419  !! integers that is safe from overflows during global
420  !! sums. This will be larger than the compile-time
421  !! precision parameter, and is used to detect overflows.
422  logical, optional, intent(inout) :: overflow !< Returns true if the conversion is being
423  !! done on a value that is too large to be represented
424  integer(kind=8), dimension(ni) :: ints
425  ! This subroutine converts a real number to an equivalent representation
426  ! using several long integers.
427 
428  real :: rs
429  character(len=80) :: mesg
430  integer(kind=8) :: ival, prec_err
431  integer :: sgn, i
432 
433  prec_err = prec ; if (present(prec_error)) prec_err = prec_error
434  ints(:) = 0_8
435  if ((r >= 1e30) .eqv. (r < 1e30)) then ; nan_error = .true. ; return ; endif
436 
437  sgn = 1 ; if (r<0.0) sgn = -1
438  rs = abs(r)
439 
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))
446  endif
447 
448  do i=1,ni
449  ival = int(rs*i_pr(i), 8)
450  rs = rs - ival*pr(i)
451  ints(i) = sgn*ival
452  enddo
453 
454 end function real_to_ints
455 
456 !> Convert the array of integers that constitute an extended-fixed-point
457 !! representation into a real number
458 function ints_to_real(ints) result(r)
459  integer(kind=8), dimension(ni), intent(in) :: ints !< The array of EFP integers
460  real :: r
461  ! This subroutine reverses the conversion in real_to_ints.
462 
463  integer :: i
464 
465  r = 0.0
466  do i=1,ni ; r = r + pr(i)*ints(i) ; enddo
467 end function ints_to_real
468 
469 !> Increment an array of integers that constitutes an extended-fixed-point
470 !! representation with a another EFP number
471 subroutine increment_ints(int_sum, int2, prec_error)
472  integer(kind=8), dimension(ni), intent(inout) :: int_sum !< The array of EFP integers being incremented
473  integer(kind=8), dimension(ni), intent(in) :: int2 !< The array of EFP integers being added
474  integer(kind=8), optional, intent(in) :: prec_error !< The PE-count dependent precision of the
475  !! integers that is safe from overflows during global
476  !! sums. This will be larger than the compile-time
477  !! precision parameter, and is used to detect overflows.
478 
479  ! This subroutine increments a number with another, both using the integer
480  ! representation in real_to_ints.
481  integer :: i
482 
483  do i=ni,2,-1
484  int_sum(i) = int_sum(i) + int2(i)
485  ! Carry the local overflow.
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
492  endif
493  enddo
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.
497  else
498  if (abs(int_sum(1)) > prec) overflow_error = .true.
499  endif
500 
501 end subroutine increment_ints
502 
503 !> Increment an EFP number with a real number without doing any carrying of
504 !! of overflows and using only minimal error checking.
505 subroutine increment_ints_faster(int_sum, r, max_mag_term)
506  integer(kind=8), dimension(ni), intent(inout) :: int_sum !< The array of EFP integers being incremented
507  real, intent(in) :: r !< The real number being added.
508  real, intent(inout) :: max_mag_term !< A running maximum magnitude of the r's.
509 
510  ! This subroutine increments a number with another, both using the integer
511  ! representation in real_to_ints, but without doing any carrying of overflow.
512  ! The entire operation is embedded in a single call for greater speed.
513  real :: rs
514  integer(kind=8) :: 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  ! Abort if the number has no EFP representation
523  if (rs > max_efp_float) then
524  overflow_error = .true.
525  return
526  endif
527 
528  do i=1,ni
529  ival = int(rs*i_pr(i), 8)
530  rs = rs - ival*pr(i)
531  int_sum(i) = int_sum(i) + sgn*ival
532  enddo
533 
534 end subroutine increment_ints_faster
535 
536 !> This subroutine handles carrying of the overflow.
537 subroutine carry_overflow(int_sum, prec_error)
538  integer(kind=8), dimension(ni), intent(inout) :: int_sum !< The array of EFP integers being
539  !! modified by carries, but without changing value.
540  integer(kind=8), intent(in) :: prec_error !< The PE-count dependent precision of the
541  !! integers that is safe from overflows during global
542  !! sums. This will be larger than the compile-time
543  !! precision parameter, and is used to detect overflows.
544 
545  ! This subroutine handles carrying of the overflow.
546  integer :: i, num_carry
547 
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
552  endif ; enddo
553  if (abs(int_sum(1)) > prec_error) then
554  overflow_error = .true.
555  endif
556 
557 end subroutine carry_overflow
558 
559 !> This subroutine carries the overflow, and then makes sure that
560 !! all integers are of the same sign as the overall value.
561 subroutine regularize_ints(int_sum)
562  integer(kind=8), dimension(ni), &
563  intent(inout) :: int_sum !< The array of integers being modified to take a
564  !! regular form with all integers of the same sign,
565  !! but without changing value.
566 
567  ! This subroutine carries the overflow, and then makes sure that
568  ! all integers are of the same sign as the overall value.
569  logical :: positive
570  integer :: i, num_carry
571 
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
576  endif ; enddo
577 
578  ! Determine the sign of the final number.
579  positive = .true.
580  do i=1,ni
581  if (abs(int_sum(i)) > 0) then
582  if (int_sum(i) < 0) positive = .false.
583  exit
584  endif
585  enddo
586 
587  if (positive) then
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
591  endif ; enddo
592  else
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
596  endif ; enddo
597  endif
598 
599 end subroutine regularize_ints
600 
601 !> Returns the status of the module's error flag
602 function query_efp_overflow_error()
605 end function query_efp_overflow_error
606 
607 !> Reset the module's error flag to false
608 subroutine reset_efp_overflow_error()
609  overflow_error = .false.
610 end subroutine reset_efp_overflow_error
611 
612 !> Add two extended-fixed-point numbers
613 function efp_plus(EFP1, EFP2)
614  type(efp_type) :: efp_plus !< The result in extended fixed point format
615  type(efp_type), intent(in) :: efp1 !< The first extended fixed point number
616  type(efp_type), intent(in) :: efp2 !< The second extended fixed point number
617 
618  efp_plus = efp1
619 
620  call increment_ints(efp_plus%v(:), efp2%v(:))
621 end function efp_plus
622 
623 !> Subract one extended-fixed-point number from another
624 function efp_minus(EFP1, EFP2)
625  type(efp_type) :: efp_minus !< The result in extended fixed point format
626  type(efp_type), intent(in) :: efp1 !< The first extended fixed point number
627  type(efp_type), intent(in) :: efp2 !< The extended fixed point number being
628  !! subtracted from the first extended fixed point number
629  integer :: i
630 
631  do i=1,ni ; efp_minus%v(i) = -1*efp2%v(i) ; enddo
632 
633  call increment_ints(efp_minus%v(:), efp1%v(:))
634 end function efp_minus
635 
636 !> Copy one extended-fixed-point number into another
637 subroutine efp_assign(EFP1, EFP2)
638  type(efp_type), intent(out) :: EFP1 !< The recipient extended fixed point number
639  type(efp_type), intent(in) :: EFP2 !< The source extended fixed point number
640  integer i
641  ! This subroutine assigns all components of the extended fixed point type
642  ! variable on the RHS (EFP2) to the components of the variable on the LHS
643  ! (EFP1).
644 
645  do i=1,ni ; efp1%v(i) = efp2%v(i) ; enddo
646 end subroutine efp_assign
647 
648 !> Return the real number that an extended-fixed-point number corresponds with
649 function efp_to_real(EFP1)
650  type(efp_type), intent(inout) :: efp1 !< The extended fixed point number being converted
651  real :: efp_to_real
652 
653  call regularize_ints(efp1%v)
654  efp_to_real = ints_to_real(efp1%v)
655 end function efp_to_real
656 
657 !> Take the difference between two extended-fixed-point numbers (EFP1 - EFP2)
658 !! and return the result as a real number
659 function efp_real_diff(EFP1, EFP2)
660  type(efp_type), intent(in) :: efp1 !< The first extended fixed point number
661  type(efp_type), intent(in) :: efp2 !< The extended fixed point number being
662  !! subtracted from the first extended fixed point number
663  real :: efp_real_diff !< The real result
664 
665  type(efp_type) :: efp_diff
666 
667  efp_diff = efp1 - efp2
668  efp_real_diff = efp_to_real(efp_diff)
669 
670 end function efp_real_diff
671 
672 !> Return the extended-fixed-point number that a real number corresponds with
673 function real_to_efp(val, overflow)
674  real, intent(in) :: val !< The real number being converted
675  logical, optional, intent(inout) :: overflow !< Returns true if the conversion is being
676  !! done on a value that is too large to be represented
677  type(efp_type) :: real_to_efp
678 
679  logical :: over
680  character(len=80) :: mesg
681 
682  if (present(overflow)) then
683  real_to_efp%v(:) = real_to_ints(val, overflow=overflow)
684  else
685  over = .false.
686  real_to_efp%v(:) = real_to_ints(val, overflow=over)
687  if (over) then
688  write(mesg, '(ES13.5)') val
689  call mom_error(fatal,"Overflow in real_to_EFP conversion of "//trim(mesg))
690  endif
691  endif
692 
693 end function real_to_efp
694 
695 !> This subroutine does a sum across PEs of a list of EFP variables,
696 !! returning the sums in place, with all overflows carried.
697 subroutine efp_list_sum_across_pes(EFPs, nval, errors)
698  type(efp_type), dimension(:), &
699  intent(inout) :: efps !< The list of extended fixed point numbers
700  !! being summed across PEs.
701  integer, intent(in) :: nval !< The number of values being summed.
702  logical, dimension(:), &
703  optional, intent(out) :: errors !< A list of error flags for each sum
704 
705  ! This subroutine does a sum across PEs of a list of EFP variables,
706  ! returning the sums in place, with all overflows carried.
707 
708  integer(kind=8), dimension(ni,nval) :: ints
709  integer(kind=8) :: prec_error
710  logical :: error_found
711  character(len=256) :: mesg
712  integer :: i, n
713 
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.")
717 
718  prec_error = (2_8**62 + (2_8**62 - 1)) / num_pes()
719  ! overflow_error is an overflow error flag for the whole module.
720  overflow_error = .false. ; error_found = .false.
721 
722  do i=1,nval ; do n=1,ni ; ints(n,i) = efps(i)%v(n) ; enddo ; enddo
723 
724  call sum_across_pes(ints(:,:), ni*nval)
725 
726  if (present(errors)) errors(:) = .false.
727  do i=1,nval
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)
736  endif
737  error_found = error_found .or. overflow_error
738  enddo
739  if (error_found .and. .not.(present(errors))) then
740  call mom_error(fatal, "Overflow in EFP_list_sum_across_PEs.")
741  endif
742 
743 end subroutine efp_list_sum_across_pes
744 
745 !> This subroutine carries out all of the calls required to close out the infrastructure cleanly.
746 !! This should only be called in ocean-only runs, as the coupler takes care of this in coupled runs.
747 subroutine mom_infra_end
748  call print_memuse_stats( 'Memory HiWaterMark', always=.true. )
749  call fms_end
750 end subroutine mom_infra_end
751 
752 end module mom_coms
mom_coms::prec
integer(kind=8), parameter prec
The precision of each integer.
Definition: MOM_coms.F90:28
mom_coms::max_efp_float
real, parameter max_efp_float
The largest float with an EFP representation. NOTE: Only the first bin can exceed precision,...
Definition: MOM_coms.F90:44
mom_coms::reproducing_sum_3d
real function reproducing_sum_3d(array, isr, ier, jsr, jer, sums, EFP_sum, err)
This subroutine uses a conversion to an integer representation of real numbers to give an order-invar...
Definition: MOM_coms.F90:245
mom_coms::max_count_prec
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: MOM_coms.F90:31
mom_coms::mom_infra_end
subroutine, public mom_infra_end
This subroutine carries out all of the calls required to close out the infrastructure cleanly....
Definition: MOM_coms.F90:748
mom_coms::efp_type
The Extended Fixed Point (EFP) type provides a public interface for doing sums and taking differences...
Definition: MOM_coms.F90:64
mom_coms::efp_plus
type(efp_type) function, public efp_plus(EFP1, EFP2)
Add two extended-fixed-point numbers.
Definition: MOM_coms.F90:614
mom_coms::overflow_error
logical overflow_error
This becomes true if an overflow is encountered.
Definition: MOM_coms.F90:49
mom_error_handler::mom_mesg
subroutine, public mom_mesg(message, verb, all_print)
This provides a convenient interface for writing an informative comment.
Definition: MOM_error_handler.F90:53
mom_coms::reproducing_sum_2d
real function reproducing_sum_2d(array, isr, ier, jsr, jer, EFP_sum, reproducing, overflow_check, err)
This subroutine uses a conversion to an integer representation of real numbers to give an order-invar...
Definition: MOM_coms.F90:83
mom_coms::query_efp_overflow_error
logical function, public query_efp_overflow_error()
Returns the status of the module's error flag.
Definition: MOM_coms.F90:603
mom_coms::efp_to_real
real function, public efp_to_real(EFP1)
Return the real number that an extended-fixed-point number corresponds with.
Definition: MOM_coms.F90:650
mom_coms::i_prec
real, parameter i_prec
The inverse of prec.
Definition: MOM_coms.F90:30
mom_coms::real_to_efp
type(efp_type) function, public real_to_efp(val, overflow)
Return the extended-fixed-point number that a real number corresponds with.
Definition: MOM_coms.F90:674
mom_coms::increment_ints
subroutine increment_ints(int_sum, int2, prec_error)
Increment an array of integers that constitutes an extended-fixed-point representation with a another...
Definition: MOM_coms.F90:472
mom_coms::carry_overflow
subroutine carry_overflow(int_sum, prec_error)
This subroutine handles carrying of the overflow.
Definition: MOM_coms.F90:538
mom_coms::regularize_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...
Definition: MOM_coms.F90:562
mom_coms
Interfaces to non-domain-oriented communication subroutines, including the MOM6 reproducing sums faci...
Definition: MOM_coms.F90:3
mom_coms::nan_error
logical nan_error
This becomes true if a NaN is encountered.
Definition: MOM_coms.F90:50
mom_coms::debug
logical debug
Making this true enables debugging output.
Definition: MOM_coms.F90:51
mom_coms::real_to_ints
integer(kind=8) function, dimension(ni) real_to_ints(r, prec_error, overflow)
Convert a real number into the array of integers constitute its extended-fixed-point representation.
Definition: MOM_coms.F90:417
mom_coms::pr
real, dimension(ni), parameter pr
An array of the real precision of each of the integers.
Definition: MOM_coms.F90:38
mom_coms::efp_real_diff
real function, public efp_real_diff(EFP1, EFP2)
Take the difference between two extended-fixed-point numbers (EFP1 - EFP2) and return the result as a...
Definition: MOM_coms.F90:660
mom_coms::i_pr
real, dimension(ni), parameter i_pr
An array of the inverse of the real precision of each of the integers.
Definition: MOM_coms.F90:41
mom_coms::increment_ints_faster
subroutine increment_ints_faster(int_sum, r, max_mag_term)
Increment an EFP number with a real number without doing any carrying of of overflows and using only ...
Definition: MOM_coms.F90:506
mom_coms::ints_to_real
real function ints_to_real(ints)
Convert the array of integers that constitute an extended-fixed-point representation into a real numb...
Definition: MOM_coms.F90:459
mom_error_handler::mom_error
subroutine, public mom_error(level, message, all_print)
This provides a convenient interface for writing an mpp_error message with run-time filter based on a...
Definition: MOM_error_handler.F90:72
mom_coms::ni
integer, parameter ni
The number of long integers to use to represent a real number.
Definition: MOM_coms.F90:36
mom_coms::reproducing_sum
Find an accurate and order-invariant sum of distributed 2d or 3d fields.
Definition: MOM_coms.F90:54
mom_coms::efp_minus
type(efp_type) function, public efp_minus(EFP1, EFP2)
Subract one extended-fixed-point number from another.
Definition: MOM_coms.F90:625
mom_coms::reset_efp_overflow_error
subroutine, public reset_efp_overflow_error()
Reset the module's error flag to false.
Definition: MOM_coms.F90:609
mom_error_handler
Routines for error handling and I/O management.
Definition: MOM_error_handler.F90:2
mom_coms::efp_list_sum_across_pes
subroutine, public 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: MOM_coms.F90:698
mom_coms::r_prec
real, parameter r_prec
A real version of prec.
Definition: MOM_coms.F90:29