MOM_debugging.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> Provides checksumming functions for debugging
6!!
7!! This module contains subroutines that perform various error checking and
8!! debugging functions for MOM6. This routine is similar to it counterpart in
9!! the SIS2 code, except for the use of the ocean_grid_type and by keeping them
10!! separate we retain the ability to set up MOM6 and SIS2 debugging separately.
11module mom_debugging
12
13use mom_checksums, only : hchksum, bchksum, qchksum, uvchksum, hchksum_pair
14use mom_checksums, only : is_nan, chksum, mom_checksums_init
15use mom_coms, only : pe_here, root_pe, num_pes
16use mom_coms, only : min_across_pes, max_across_pes, reproducing_sum
17use mom_domains, only : pass_vector, pass_var, pe_here
18use mom_domains, only : bgrid_ne, agrid, to_all, scalar_pair
19use mom_error_handler, only : mom_error, fatal, warning, is_root_pe
20use mom_file_parser, only : log_version, param_file_type, get_param
21use mom_grid, only : ocean_grid_type
23use mom_io, only : stdout
25
26implicit none ; private
27
28public :: check_redundant_c, check_redundant_b, check_redundant_t, check_redundant
29public :: vec_chksum, vec_chksum_c, vec_chksum_b, vec_chksum_a
33
34! These interfaces come from MOM_checksums.
35public :: hchksum, bchksum, qchksum, is_nan, chksum, uvchksum, hchksum_pair
36
37!> Check for consistency between the duplicated points of a C-grid vector
38interface check_redundant
40end interface check_redundant
41!> Check for consistency between the duplicated points of a C-grid vector
42interface check_redundant_c
44end interface check_redundant_c
45!> Check for consistency between the duplicated points of a B-grid vector or scalar
46interface check_redundant_b
49end interface check_redundant_b
50!> Check for consistency between the duplicated points of an A-grid vector or scalar
51interface check_redundant_t
54end interface check_redundant_t
55
56!> Do checksums on the components of a C-grid vector
57interface vec_chksum
58 module procedure chksum_vec_c3d, chksum_vec_c2d
59end interface vec_chksum
60!> Do checksums on the components of a C-grid vector
61interface vec_chksum_c
62 module procedure chksum_vec_c3d, chksum_vec_c2d
63end interface vec_chksum_c
64!> Do checksums on the components of a B-grid vector
65interface vec_chksum_b
66 module procedure chksum_vec_b3d, chksum_vec_b2d
67end interface vec_chksum_b
68!> Do checksums on the components of an A-grid vector
69interface vec_chksum_a
70 module procedure chksum_vec_a3d, chksum_vec_a2d
71end interface vec_chksum_a
72
73! Note: these parameters are module data but ONLY used when debugging and
74! so can violate the thread-safe requirement of no module/global data.
75integer :: max_redundant_prints = 100 !< Maximum number of times to write redundant messages
76integer :: redundant_prints(3) = 0 !< Counters for controlling redundant printing
77logical :: debug = .false. !< Write out verbose debugging data
78logical :: debug_chksums = .true. !< Perform checksums on arrays
79logical :: debug_redundant = .true. !< Check redundant values on PE boundaries
80
81contains
82
83!> MOM_debugging_init initializes the MOM_debugging module, and sets
84!! the parameters that control which checks are active for MOM6.
85subroutine mom_debugging_init(param_file)
86 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
87 ! This include declares and sets the variable "version".
88# include "version_variable.h"
89 character(len=40) :: mdl = "MOM_debugging" ! This module's name.
90
91 call log_version(param_file, mdl, version, debugging=.true.)
92 call get_param(param_file, mdl, "DEBUG", debug, &
93 "If true, write out verbose debugging data.", &
94 default=.false., debuggingparam=.true.)
95 call get_param(param_file, mdl, "DEBUG_CHKSUMS", debug_chksums, &
96 "If true, checksums are performed on arrays in the "//&
97 "various vec_chksum routines.", default=debug, &
98 debuggingparam=.true.)
99 call get_param(param_file, mdl, "DEBUG_REDUNDANT", debug_redundant, &
100 "If true, debug redundant data points during calls to "//&
101 "the various vec_chksum routines.", default=debug, &
102 debuggingparam=.true.)
103
104 call mom_checksums_init(param_file)
105
106end subroutine mom_debugging_init
107
108!> Returns logicals indicating which debugging checks should be performed.
109subroutine query_debugging_checks(do_debug, do_chksums, do_redundant)
110 logical, optional, intent(out) :: do_debug !< True if verbose debugging is to be output
111 logical, optional, intent(out) :: do_chksums !< True if checksums are to be output
112 logical, optional, intent(out) :: do_redundant !< True if redundant points are to be checked
113
114 if (present(do_debug)) do_debug = debug
115 if (present(do_chksums)) do_chksums = debug_chksums
116 if (present(do_redundant)) do_redundant = debug_redundant
117
118end subroutine query_debugging_checks
119
120!> Check for consistency between the duplicated points of a 3-D C-grid vector
121subroutine check_redundant_vc3d(mesg, u_comp, v_comp, G, is, ie, js, je, &
122 direction, unscale)
123 character(len=*), intent(in) :: mesg !< An identifying message
124 type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
125 real, dimension(G%IsdB:,G%jsd:,:), intent(in) :: u_comp !< The u-component of the vector to be
126 !! checked for consistency in arbitrary,
127 !! possibly rescaled units [A ~> a]
128 real, dimension(G%isd:,G%JsdB:,:), intent(in) :: v_comp !< The u-component of the vector to be
129 !! checked for consistency in arbitrary,
130 !! possibly rescaled units [A ~> a]
131 integer, optional, intent(in) :: is !< The starting i-index to check
132 integer, optional, intent(in) :: ie !< The ending i-index to check
133 integer, optional, intent(in) :: js !< The starting j-index to check
134 integer, optional, intent(in) :: je !< The ending j-index to check
135 integer, optional, intent(in) :: direction !< the direction flag to be
136 !! passed to pass_vector
137 real, optional, intent(in) :: unscale !< A factor that undoes the scaling for the
138 !! arrays to give consistent output [a A-1 ~> 1]
139
140 ! Local variables
141 character(len=24) :: mesg_k
142 integer :: k
143
144 do k=1,size(u_comp,3)
145 write(mesg_k,'(" Layer ",i0," ")') k
146 call check_redundant_vc2d(trim(mesg)//trim(mesg_k), u_comp(:,:,k), &
147 v_comp(:,:,k), g, is, ie, js, je, direction, unscale)
148 enddo
149end subroutine check_redundant_vc3d
150
151!> Check for consistency between the duplicated points of a 2-D C-grid vector
152subroutine check_redundant_vc2d(mesg, u_comp, v_comp, G, is, ie, js, je, &
153 direction, unscale)
154 character(len=*), intent(in) :: mesg !< An identifying message
155 type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
156 real, dimension(G%IsdB:,G%jsd:), intent(in) :: u_comp !< The u-component of the vector to be
157 !! checked for consistency in arbitrary,
158 !! possibly rescaled units [A ~> a]
159 real, dimension(G%isd:,G%JsdB:), intent(in) :: v_comp !< The u-component of the vector to be
160 !! checked for consistency in arbitrary,
161 !! possibly rescaled units [A ~> a]
162 integer, optional, intent(in) :: is !< The starting i-index to check
163 integer, optional, intent(in) :: ie !< The ending i-index to check
164 integer, optional, intent(in) :: js !< The starting j-index to check
165 integer, optional, intent(in) :: je !< The ending j-index to check
166 integer, optional, intent(in) :: direction !< the direction flag to be
167 !! passed to pass_vector
168 real, optional, intent(in) :: unscale !< A factor that undoes the scaling for the
169 !! arrays to give consistent output [a A-1 ~> 1]
170 ! Local variables
171 ! In the following comments, [A] is used to indicate the arbitrary, possibly rescaled units
172 ! of the input vector while [a] indicates the unscaled (e.g., mks) units to used for output.
173 real :: u_nonsym(G%isd:G%ied,G%jsd:G%jed) ! A nonsymmetric version of u_comp [A ~> a]
174 real :: v_nonsym(G%isd:G%ied,G%jsd:G%jed) ! A nonsymmetric version of v_comp [A ~> a]
175 real :: u_resym(G%IsdB:G%IedB,G%jsd:G%jed) ! A reconstructed symmetric version of u_comp [A ~> a]
176 real :: v_resym(G%isd:G%ied,G%JsdB:G%JedB) ! A reconstructed symmetric version of v_comp [A ~> a]
177 real :: sc ! A factor that undoes the scaling for the arrays to give consistent output [a A-1 ~> 1]
178 character(len=128) :: mesg2
179 integer :: i, j, is_ch, ie_ch, js_ch, je_ch
180 integer :: Isq, Ieq, Jsq, Jeq, isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
181
182 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
183 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
184 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
185
186 if (.not.(present(is) .or. present(ie) .or. present(js) .or. present(je))) then
187 ! This only works with symmetric memory, so otherwise return.
188 if ((isd == isdb) .and. (jsd == jsdb)) return
189 endif
190
191 sc = 1.0 ; if (present(unscale)) sc = unscale
192
193 do i=isd,ied ; do j=jsd,jed
194 u_nonsym(i,j) = u_comp(i,j) ; v_nonsym(i,j) = v_comp(i,j)
195 enddo ; enddo
196
197 if (.not.associated(g%Domain_aux)) call mom_error(fatal, &
198 " check_redundant called with a non-associated auxiliary domain the grid type.")
199 call pass_vector(u_nonsym, v_nonsym, g%Domain_aux, direction)
200
201 do i=isdb,iedb ; do j=jsd,jed ; u_resym(i,j) = u_comp(i,j) ; enddo ; enddo
202 do i=isd,ied ; do j=jsdb,jedb ; v_resym(i,j) = v_comp(i,j) ; enddo ; enddo
203 do i=isd,ied ; do j=jsd,jed
204 u_resym(i,j) = u_nonsym(i,j) ; v_resym(i,j) = v_nonsym(i,j)
205 enddo ; enddo
206 call pass_vector(u_resym, v_resym, g%Domain, direction)
207
208 is_ch = isq ; ie_ch = ieq ; js_ch = jsq ; je_ch = jeq
209 if (present(is)) is_ch = is ; if (present(ie)) ie_ch = ie
210 if (present(js)) js_ch = js ; if (present(js)) je_ch = je
211
212 do i=is_ch,ie_ch ; do j=js_ch+1,je_ch
213 if (u_resym(i,j) /= u_comp(i,j) .and. &
214 redundant_prints(3) < max_redundant_prints) then
215 write(mesg2,'(" redundant u-components",2(1pe12.4)," differ by ", &
216 & 1pe12.4," at i,j = ",I0,",",I0," on pe ",I0)') &
217 sc*u_comp(i,j), sc*u_resym(i,j), sc*(u_comp(i,j)-u_resym(i,j)), i, j, pe_here()
218 write(0,'(A130)') trim(mesg)//trim(mesg2)
219 redundant_prints(3) = redundant_prints(3) + 1
220 endif
221 enddo ; enddo
222 do i=is_ch+1,ie_ch ; do j=js_ch,je_ch
223 if (v_resym(i,j) /= v_comp(i,j) .and. &
224 redundant_prints(3) < max_redundant_prints) then
225 write(mesg2,'(" redundant v-comps",2(1pe12.4)," differ by ", &
226 & 1pe12.4," at i,j = ",I0,",",I0," x,y = ",2(1pe12.4)," on pe ",I0)') &
227 sc*v_comp(i,j), sc*v_resym(i,j), sc*(v_comp(i,j)-v_resym(i,j)), i, j, &
228 g%geoLonBu(i,j), g%geoLatBu(i,j), pe_here()
229 write(0,'(A155)') trim(mesg)//trim(mesg2)
230 redundant_prints(3) = redundant_prints(3) + 1
231 endif
232 enddo ; enddo
233
234end subroutine check_redundant_vc2d
235
236!> Check for consistency between the duplicated points of a 3-D scalar at corner points
237subroutine check_redundant_sb3d(mesg, array, G, is, ie, js, je, unscale)
238 character(len=*), intent(in) :: mesg !< An identifying message
239 type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
240 real, dimension(G%IsdB:,G%JsdB:,:), intent(in) :: array !< The array to be checked for consistency in
241 !! arbitrary, possibly rescaled units [A ~> a]
242 integer, optional, intent(in) :: is !< The starting i-index to check
243 integer, optional, intent(in) :: ie !< The ending i-index to check
244 integer, optional, intent(in) :: js !< The starting j-index to check
245 integer, optional, intent(in) :: je !< The ending j-index to check
246 real, optional, intent(in) :: unscale !< A factor that undoes the scaling for the
247 !! arrays to give consistent output [a A-1 ~> 1]
248
249 ! Local variables
250 character(len=24) :: mesg_k
251 integer :: k
252
253 do k=1,size(array,3)
254 write(mesg_k,'(" Layer ",i0," ")') k
255 call check_redundant_sb2d(trim(mesg)//trim(mesg_k), array(:,:,k), &
256 g, is, ie, js, je, unscale)
257 enddo
258end subroutine check_redundant_sb3d
259
260!> Check for consistency between the duplicated points of a 2-D scalar at corner points
261subroutine check_redundant_sb2d(mesg, array, G, is, ie, js, je, unscale)
262 character(len=*), intent(in) :: mesg !< An identifying message
263 type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
264 real, dimension(G%IsdB:,G%JsdB:), intent(in) :: array !< The array to be checked for consistency in
265 !! arbitrary, possibly rescaled units [A ~> a]
266 integer, optional, intent(in) :: is !< The starting i-index to check
267 integer, optional, intent(in) :: ie !< The ending i-index to check
268 integer, optional, intent(in) :: js !< The starting j-index to check
269 integer, optional, intent(in) :: je !< The ending j-index to check
270 real, optional, intent(in) :: unscale !< A factor that undoes the scaling for the
271 !! arrays to give consistent output [a A-1 ~> 1]
272 ! Local variables
273 ! In the following comments, [A] is used to indicate the arbitrary, possibly rescaled units
274 ! of the input array while [a] indicates the unscaled (e.g., mks) units to used for output.
275 real :: a_nonsym(G%isd:G%ied,G%jsd:G%jed) ! A nonsymmetric version of array [A ~> a]
276 real :: a_resym(G%IsdB:G%IedB,G%JsdB:G%JedB) ! A reconstructed symmetric version of array [A ~> a]
277 real :: sc ! A factor that undoes the scaling for the arrays to give consistent output [a A-1 ~> 1]
278 character(len=128) :: mesg2
279 integer :: i, j, is_ch, ie_ch, js_ch, je_ch
280 integer :: Isq, Ieq, Jsq, Jeq, isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
281
282 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
283 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
284 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
285
286 if (.not.(present(is) .or. present(ie) .or. present(js) .or. present(je))) then
287 ! This only works with symmetric memory, so otherwise return.
288 if ((isd == isdb) .and. (jsd == jsdb)) return
289 endif
290
291 sc = 1.0 ; if (present(unscale)) sc = unscale
292
293 do i=isd,ied ; do j=jsd,jed
294 a_nonsym(i,j) = array(i,j)
295 enddo ; enddo
296
297 if (.not.associated(g%Domain_aux)) call mom_error(fatal, &
298 " check_redundant called with a non-associated auxiliary domain the grid type.")
299 call pass_vector(a_nonsym, a_nonsym, g%Domain_aux, &
300 direction=to_all+scalar_pair, stagger=bgrid_ne)
301
302 do i=isdb,iedb ; do j=jsdb,jedb ; a_resym(i,j) = array(i,j) ; enddo ; enddo
303 do i=isd,ied ; do j=jsd,jed
304 a_resym(i,j) = a_nonsym(i,j)
305 enddo ; enddo
306 call pass_vector(a_resym, a_resym, g%Domain, direction=to_all+scalar_pair, &
307 stagger=bgrid_ne)
308
309 is_ch = isq ; ie_ch = ieq ; js_ch = jsq ; je_ch = jeq
310 if (present(is)) is_ch = is ; if (present(ie)) ie_ch = ie
311 if (present(js)) js_ch = js ; if (present(js)) je_ch = je
312
313 do i=is_ch,ie_ch ; do j=js_ch,je_ch
314 if (a_resym(i,j) /= array(i,j) .and. &
315 redundant_prints(2) < max_redundant_prints) then
316 write(mesg2,'(" Redundant points",2(1pe12.4)," differ by ", &
317 & 1pe12.4," at i,j = ",I0,",",I0," on pe ",I0)') &
318 sc*array(i,j), sc*a_resym(i,j), sc*(array(i,j)-a_resym(i,j)), i, j, pe_here()
319 write(0,'(A130)') trim(mesg)//trim(mesg2)
320 redundant_prints(2) = redundant_prints(2) + 1
321 endif
322 enddo ; enddo
323
324end subroutine check_redundant_sb2d
325
326!> Check for consistency between the duplicated points of a 3-D B-grid vector
327subroutine check_redundant_vb3d(mesg, u_comp, v_comp, G, is, ie, js, je, &
328 direction, unscale)
329 character(len=*), intent(in) :: mesg !< An identifying message
330 type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
331 real, dimension(G%IsdB:,G%JsdB:,:), intent(in) :: u_comp !< The u-component of the vector to be
332 !! checked for consistency in arbitrary,
333 !! possibly rescaled units [A ~> a]
334 real, dimension(G%IsdB:,G%JsdB:,:), intent(in) :: v_comp !< The v-component of the vector to be
335 !! checked for consistency in arbitrary,
336 !! possibly rescaled units [A ~> a]
337 integer, optional, intent(in) :: is !< The starting i-index to check
338 integer, optional, intent(in) :: ie !< The ending i-index to check
339 integer, optional, intent(in) :: js !< The starting j-index to check
340 integer, optional, intent(in) :: je !< The ending j-index to check
341 integer, optional, intent(in) :: direction !< the direction flag to be
342 !! passed to pass_vector
343 real, optional, intent(in) :: unscale !< A factor that undoes the scaling for the
344 !! arrays to give consistent output [a A-1 ~> 1]
345 ! Local variables
346 character(len=24) :: mesg_k
347 integer :: k
348
349 do k=1,size(u_comp,3)
350 write(mesg_k,'(" Layer ",i0," ")') k
351 call check_redundant_vb2d(trim(mesg)//trim(mesg_k), u_comp(:,:,k), &
352 v_comp(:,:,k), g, is, ie, js, je, direction, unscale)
353 enddo
354end subroutine check_redundant_vb3d
355
356!> Check for consistency between the duplicated points of a 2-D B-grid vector
357subroutine check_redundant_vb2d(mesg, u_comp, v_comp, G, is, ie, js, je, &
358 direction, unscale)
359 character(len=*), intent(in) :: mesg !< An identifying message
360 type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
361 real, dimension(G%IsdB:,G%JsdB:), intent(in) :: u_comp !< The u-component of the vector to be
362 !! checked for consistency in arbitrary,
363 !! possibly rescaled units [A ~> a]
364 real, dimension(G%IsdB:,G%JsdB:), intent(in) :: v_comp !< The v-component of the vector to be
365 !! checked for consistency in arbitrary,
366 !! possibly rescaled units [A ~> a]
367 integer, optional, intent(in) :: is !< The starting i-index to check
368 integer, optional, intent(in) :: ie !< The ending i-index to check
369 integer, optional, intent(in) :: js !< The starting j-index to check
370 integer, optional, intent(in) :: je !< The ending j-index to check
371 integer, optional, intent(in) :: direction !< the direction flag to be
372 !! passed to pass_vector
373 real, optional, intent(in) :: unscale !< A factor that undoes the scaling for the
374 !! arrays to give consistent output [a A-1 ~> 1]
375 ! Local variables
376 ! In the following comments, [A] is used to indicate the arbitrary, possibly rescaled units
377 ! of the input vector while [a] indicates the unscaled (e.g., mks) units to used for output.
378 real :: u_nonsym(G%isd:G%ied,G%jsd:G%jed) ! A nonsymmetric version of u_comp [A ~> a]
379 real :: v_nonsym(G%isd:G%ied,G%jsd:G%jed) ! A nonsymmetric version of v_comp [A ~> a]
380 real :: u_resym(G%IsdB:G%IedB,G%JsdB:G%JedB) ! A reconstructed symmetric version of u_comp [A ~> a]
381 real :: v_resym(G%IsdB:G%IedB,G%JsdB:G%JedB) ! A reconstructed symmetric version of v_comp [A ~> a]
382 real :: sc ! A factor that undoes the scaling for the arrays to give consistent output [a A-1 ~> 1]
383 character(len=128) :: mesg2
384 integer :: i, j, is_ch, ie_ch, js_ch, je_ch
385 integer :: Isq, Ieq, Jsq, Jeq, isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
386
387 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
388 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
389 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
390
391 if (.not.(present(is) .or. present(ie) .or. present(js) .or. present(je))) then
392 ! This only works with symmetric memory, so otherwise return.
393 if ((isd == isdb) .and. (jsd == jsdb)) return
394 endif
395
396 sc = 1.0 ; if (present(unscale)) sc = unscale
397
398 do i=isd,ied ; do j=jsd,jed
399 u_nonsym(i,j) = u_comp(i,j) ; v_nonsym(i,j) = v_comp(i,j)
400 enddo ; enddo
401
402 if (.not.associated(g%Domain_aux)) call mom_error(fatal, &
403 " check_redundant called with a non-associated auxiliary domain the grid type.")
404 call pass_vector(u_nonsym, v_nonsym, g%Domain_aux, direction, stagger=bgrid_ne)
405
406 do i=isdb,iedb ; do j=jsdb,jedb
407 u_resym(i,j) = u_comp(i,j) ; v_resym(i,j) = v_comp(i,j)
408 enddo ; enddo
409 do i=isd,ied ; do j=jsd,jed
410 u_resym(i,j) = u_nonsym(i,j) ; v_resym(i,j) = v_nonsym(i,j)
411 enddo ; enddo
412 call pass_vector(u_resym, v_resym, g%Domain, direction, stagger=bgrid_ne)
413
414 is_ch = isq ; ie_ch = ieq ; js_ch = jsq ; je_ch = jeq
415 if (present(is)) is_ch = is ; if (present(ie)) ie_ch = ie
416 if (present(js)) js_ch = js ; if (present(js)) je_ch = je
417
418 do i=is_ch,ie_ch ; do j=js_ch,je_ch
419 if (u_resym(i,j) /= u_comp(i,j) .and. &
420 redundant_prints(2) < max_redundant_prints) then
421 write(mesg2,'(" redundant u-components",2(1pe12.4)," differ by ", &
422 & 1pe12.4," at i,j = ",I0,",",I0," on pe ",I0)') &
423 sc*u_comp(i,j), sc*u_resym(i,j), sc*(u_comp(i,j)-u_resym(i,j)), i, j, pe_here()
424 write(0,'(A130)') trim(mesg)//trim(mesg2)
425 redundant_prints(2) = redundant_prints(2) + 1
426 endif
427 enddo ; enddo
428 do i=is_ch,ie_ch ; do j=js_ch,je_ch
429 if (v_resym(i,j) /= v_comp(i,j) .and. &
430 redundant_prints(2) < max_redundant_prints) then
431 write(mesg2,'(" redundant v-comps",2(1pe12.4)," differ by ", &
432 & 1pe12.4," at i,j = ",I0,",",I0," x,y = ",2(1pe12.4)," on pe ",I0)') &
433 sc*v_comp(i,j), sc*v_resym(i,j), sc*(v_comp(i,j)-v_resym(i,j)), i, j, &
434 g%geoLonBu(i,j), g%geoLatBu(i,j), pe_here()
435 write(0,'(A155)') trim(mesg)//trim(mesg2)
436 redundant_prints(2) = redundant_prints(2) + 1
437 endif
438 enddo ; enddo
439
440end subroutine check_redundant_vb2d
441
442!> Check for consistency between the duplicated points of a 3-D scalar at tracer points
443subroutine check_redundant_st3d(mesg, array, G, is, ie, js, je, unscale)
444 character(len=*), intent(in) :: mesg !< An identifying message
445 type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
446 real, dimension(G%isd:,G%jsd:,:), intent(in) :: array !< The array to be checked for consistency in
447 !! arbitrary, possibly rescaled units [A ~> a]
448 integer, optional, intent(in) :: is !< The starting i-index to check
449 integer, optional, intent(in) :: ie !< The ending i-index to check
450 integer, optional, intent(in) :: js !< The starting j-index to check
451 integer, optional, intent(in) :: je !< The ending j-index to check
452 real, optional, intent(in) :: unscale !< A factor that undoes the scaling for the
453 !! arrays to give consistent output [a A-1 ~> 1]
454 ! Local variables
455 character(len=24) :: mesg_k
456 integer :: k
457
458 do k=1,size(array,3)
459 write(mesg_k,'(" Layer ",i0," ")') k
460 call check_redundant_st2d(trim(mesg)//trim(mesg_k), array(:,:,k), &
461 g, is, ie, js, je, unscale)
462 enddo
463end subroutine check_redundant_st3d
464
465
466!> Check for consistency between the duplicated points of a 2-D scalar at tracer points
467subroutine check_redundant_st2d(mesg, array, G, is, ie, js, je, unscale)
468 character(len=*), intent(in) :: mesg !< An identifying message
469 type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
470 real, dimension(G%isd:,G%jsd:), intent(in) :: array !< The array to be checked for consistency in
471 !! arbitrary, possibly rescaled units [A ~> a]
472 integer, optional, intent(in) :: is !< The starting i-index to check
473 integer, optional, intent(in) :: ie !< The ending i-index to check
474 integer, optional, intent(in) :: js !< The starting j-index to check
475 integer, optional, intent(in) :: je !< The ending j-index to check
476 real, optional, intent(in) :: unscale !< A factor that undoes the scaling for the
477 !! arrays to give consistent output [a A-1 ~> 1]
478 ! Local variables
479 ! In the following comments, [A] is used to indicate the arbitrary, possibly rescaled units
480 ! of the input array while [a] indicates the unscaled (e.g., mks) units to used for output.
481 real :: a_nonsym(G%isd:G%ied,G%jsd:G%jed) ! A version of array with halo points updated by message passing [A ~> a]
482 real :: sc ! A factor that undoes the scaling for the arrays to give consistent output [a A-1 ~> 1]
483 character(len=128) :: mesg2
484
485 integer :: i, j, is_ch, ie_ch, js_ch, je_ch
486 integer :: isd, ied, jsd, jed
487 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
488
489 is_ch = g%isc ; ie_ch = g%iec ; js_ch = g%jsc ; je_ch = g%jec
490 if (present(is)) is_ch = is ; if (present(ie)) ie_ch = ie
491 if (present(js)) js_ch = js ; if (present(js)) je_ch = je
492
493 sc = 1.0 ; if (present(unscale)) sc = unscale
494
495 ! This only works on points outside of the standard computational domain.
496 if ((is_ch == g%isc) .and. (ie_ch == g%iec) .and. &
497 (js_ch == g%jsc) .and. (je_ch == g%jec)) return
498
499 do i=isd,ied ; do j=jsd,jed
500 a_nonsym(i,j) = array(i,j)
501 enddo ; enddo
502
503 call pass_var(a_nonsym, g%Domain)
504
505 do i=is_ch,ie_ch ; do j=js_ch,je_ch
506 if (a_nonsym(i,j) /= array(i,j) .and. &
507 redundant_prints(1) < max_redundant_prints) then
508 write(mesg2,'(" Redundant points",2(1pe12.4)," differ by ", &
509 & 1pe12.4," at i,j = ",I0,",",I0," on pe ",I0)') &
510 sc*array(i,j), sc*a_nonsym(i,j), sc*(array(i,j)-a_nonsym(i,j)), i, j, pe_here()
511 write(0,'(A130)') trim(mesg)//trim(mesg2)
512 redundant_prints(1) = redundant_prints(1) + 1
513 endif
514 enddo ; enddo
515
516end subroutine check_redundant_st2d
517
518!> Check for consistency between the duplicated points of a 3-D A-grid vector
519subroutine check_redundant_vt3d(mesg, u_comp, v_comp, G, is, ie, js, je, &
520 direction, unscale)
521 character(len=*), intent(in) :: mesg !< An identifying message
522 type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
523 real, dimension(G%isd:,G%jsd:,:), intent(in) :: u_comp !< The u-component of the vector to be
524 !! checked for consistency in arbitrary,
525 !! possibly rescaled units [A ~> a]
526 real, dimension(G%isd:,G%jsd:,:), intent(in) :: v_comp !< The v-component of the vector to be
527 !! checked for consistency in arbitrary,
528 !! possibly rescaled units [A ~> a]
529 integer, optional, intent(in) :: is !< The starting i-index to check
530 integer, optional, intent(in) :: ie !< The ending i-index to check
531 integer, optional, intent(in) :: js !< The starting j-index to check
532 integer, optional, intent(in) :: je !< The ending j-index to check
533 integer, optional, intent(in) :: direction !< the direction flag to be
534 !! passed to pass_vector
535 real, optional, intent(in) :: unscale !< A factor that undoes the scaling for the
536 !! arrays to give consistent output [a A-1 ~> 1]
537 ! Local variables
538 character(len=24) :: mesg_k
539 integer :: k
540
541 do k=1,size(u_comp,3)
542 write(mesg_k,'(" Layer ",i0," ")') k
543 call check_redundant_vt2d(trim(mesg)//trim(mesg_k), u_comp(:,:,k), &
544 v_comp(:,:,k), g, is, ie, js, je, direction, unscale)
545 enddo
546end subroutine check_redundant_vt3d
547
548!> Check for consistency between the duplicated points of a 2-D A-grid vector
549subroutine check_redundant_vt2d(mesg, u_comp, v_comp, G, is, ie, js, je, &
550 direction, unscale)
551 character(len=*), intent(in) :: mesg !< An identifying message
552 type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
553 real, dimension(G%isd:,G%jsd:), intent(in) :: u_comp !< The u-component of the vector to be
554 !! checked for consistency in arbitrary,
555 !! possibly rescaled units [A ~> a]
556 real, dimension(G%isd:,G%jsd:), intent(in) :: v_comp !< The v-component of the vector to be
557 !! checked for consistency in arbitrary,
558 !! possibly rescaled units [A ~> a]
559 integer, optional, intent(in) :: is !< The starting i-index to check
560 integer, optional, intent(in) :: ie !< The ending i-index to check
561 integer, optional, intent(in) :: js !< The starting j-index to check
562 integer, optional, intent(in) :: je !< The ending j-index to check
563 integer, optional, intent(in) :: direction !< the direction flag to be
564 !! passed to pass_vector
565 real, optional, intent(in) :: unscale !< A factor that undoes the scaling for the
566 !! arrays to give consistent output [a A-1 ~> 1]
567 ! Local variables
568 ! In the following comments, [A] is used to indicate the arbitrary, possibly rescaled units
569 ! of the input vector while [a] indicates the unscaled (e.g., mks) units to used for output.
570 real :: u_nonsym(G%isd:G%ied,G%jsd:G%jed) ! A version of u_comp with halo points updated by message passing [A ~> a]
571 real :: v_nonsym(G%isd:G%ied,G%jsd:G%jed) ! A version of v_comp with halo points updated by message passing [A ~> a]
572 real :: sc ! A factor that undoes the scaling for the arrays to give consistent output [a A-1 ~> 1]
573 character(len=128) :: mesg2
574
575 integer :: i, j, is_ch, ie_ch, js_ch, je_ch
576 integer :: Isq, Ieq, Jsq, Jeq, isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
577 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
578 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
579 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
580
581 is_ch = g%isc ; ie_ch = g%iec ; js_ch = g%jsc ; je_ch = g%jec
582 if (present(is)) is_ch = is ; if (present(ie)) ie_ch = ie
583 if (present(js)) js_ch = js ; if (present(js)) je_ch = je
584
585 sc = 1.0 ; if (present(unscale)) sc = unscale
586
587 ! This only works on points outside of the standard computational domain.
588 if ((is_ch == g%isc) .and. (ie_ch == g%iec) .and. &
589 (js_ch == g%jsc) .and. (je_ch == g%jec)) return
590
591 do i=isd,ied ; do j=jsd,jed
592 u_nonsym(i,j) = u_comp(i,j) ; v_nonsym(i,j) = v_comp(i,j)
593 enddo ; enddo
594
595 call pass_vector(u_nonsym, v_nonsym, g%Domain, direction, stagger=agrid)
596
597 do i=is_ch,ie_ch ; do j=js_ch+1,je_ch
598 if (u_nonsym(i,j) /= u_comp(i,j) .and. &
599 redundant_prints(1) < max_redundant_prints) then
600 write(mesg2,'(" redundant u-components",2(1pe12.4)," differ by ", &
601 & 1pe12.4," at i,j = ",I0,",",I0," on pe ",I0)') &
602 sc*u_comp(i,j), sc*u_nonsym(i,j), sc*(u_comp(i,j)-u_nonsym(i,j)), i, j, pe_here()
603 write(0,'(A130)') trim(mesg)//trim(mesg2)
604 redundant_prints(1) = redundant_prints(1) + 1
605 endif
606 enddo ; enddo
607 do i=is_ch+1,ie_ch ; do j=js_ch,je_ch
608 if (v_nonsym(i,j) /= v_comp(i,j) .and. &
609 redundant_prints(1) < max_redundant_prints) then
610 write(mesg2,'(" redundant v-comps",2(1pe12.4)," differ by ", &
611 & 1pe12.4," at i,j = ",I0,",",I0," x,y = ",2(1pe12.4)," on pe ",I0)') &
612 sc*v_comp(i,j), sc*v_nonsym(i,j), sc*(v_comp(i,j)-v_nonsym(i,j)), i, j, &
613 g%geoLonBu(i,j), g%geoLatBu(i,j), pe_here()
614 write(0,'(A155)') trim(mesg)//trim(mesg2)
615 redundant_prints(1) = redundant_prints(1) + 1
616 endif
617 enddo ; enddo
618
619end subroutine check_redundant_vt2d
620
621
622! It appears that none of the other routines in this file are ever called.
623
624!> Do a checksum and redundant point check on a 3d C-grid vector.
625subroutine chksum_vec_c3d(mesg, u_comp, v_comp, G, halos, scalars, unscale)
626 character(len=*), intent(in) :: mesg !< An identifying message
627 type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
628 real, dimension(G%IsdB:,G%jsd:,:), intent(in) :: u_comp !< The u-component of the vector to be
629 !! checked for consistency in arbitrary,
630 !! possibly rescaled units [A ~> a]
631 real, dimension(G%isd:,G%JsdB:,:), intent(in) :: v_comp !< The v-component of the vector to be
632 !! checked for consistency in arbitrary,
633 !! possibly rescaled units [A ~> a]
634 integer, optional, intent(in) :: halos !< The width of halos to check (default 0)
635 logical, optional, intent(in) :: scalars !< If true this is a pair of
636 !! scalars that are being checked.
637 real, optional, intent(in) :: unscale !< A factor that undoes the scaling for the
638 !! arrays to give consistent output [a A-1 ~> 1]
639 ! Local variables
640 logical :: are_scalars
641 are_scalars = .false. ; if (present(scalars)) are_scalars = scalars
642
643 if (debug_chksums) then
644 call uvchksum(mesg, u_comp, v_comp, g%HI, halos, unscale=unscale)
645 endif
646 if (debug_redundant) then
647 if (are_scalars) then
648 call check_redundant_c(mesg, u_comp, v_comp, g, direction=to_all+scalar_pair, unscale=unscale)
649 else
650 call check_redundant_c(mesg, u_comp, v_comp, g, unscale=unscale)
651 endif
652 endif
653
654end subroutine chksum_vec_c3d
655
656!> Do a checksum and redundant point check on a 2d C-grid vector.
657subroutine chksum_vec_c2d(mesg, u_comp, v_comp, G, halos, scalars, unscale)
658 character(len=*), intent(in) :: mesg !< An identifying message
659 type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
660 real, dimension(G%IsdB:,G%jsd:), intent(in) :: u_comp !< The u-component of the vector to be
661 !! checked for consistency in arbitrary,
662 !! possibly rescaled units [A ~> a]
663 real, dimension(G%isd:,G%JsdB:), intent(in) :: v_comp !< The v-component of the vector to be
664 !! checked for consistency in arbitrary,
665 !! possibly rescaled units [A ~> a]
666 integer, optional, intent(in) :: halos !< The width of halos to check (default 0)
667 logical, optional, intent(in) :: scalars !< If true this is a pair of
668 !! scalars that are being checked.
669 real, optional, intent(in) :: unscale !< A factor that undoes the scaling for the
670 !! arrays to give consistent output [a A-1 ~> 1]
671 ! Local variables
672 logical :: are_scalars
673 are_scalars = .false. ; if (present(scalars)) are_scalars = scalars
674
675 if (debug_chksums) then
676 call uvchksum(mesg, u_comp, v_comp, g%HI, halos, unscale=unscale)
677 endif
678 if (debug_redundant) then
679 if (are_scalars) then
680 call check_redundant_c(mesg, u_comp, v_comp, g, direction=to_all+scalar_pair, unscale=unscale)
681 else
682 call check_redundant_c(mesg, u_comp, v_comp, g, unscale=unscale)
683 endif
684 endif
685
686end subroutine chksum_vec_c2d
687
688!> Do a checksum and redundant point check on a 3d B-grid vector.
689subroutine chksum_vec_b3d(mesg, u_comp, v_comp, G, halos, scalars, unscale)
690 character(len=*), intent(in) :: mesg !< An identifying message
691 type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
692 real, dimension(G%IsdB:,G%JsdB:,:), intent(in) :: u_comp !< The u-component of the vector to be
693 !! checked for consistency in arbitrary,
694 !! possibly rescaled units [A ~> a]
695 real, dimension(G%IsdB:,G%JsdB:,:), intent(in) :: v_comp !< The v-component of the vector to be
696 !! checked for consistency in arbitrary,
697 !! possibly rescaled units [A ~> a]
698 integer, optional, intent(in) :: halos !< The width of halos to check (default 0)
699 logical, optional, intent(in) :: scalars !< If true this is a pair of
700 !! scalars that are being checked.
701 real, optional, intent(in) :: unscale !< A factor that undoes the scaling for the
702 !! arrays to give consistent output [a A-1 ~> 1]
703 ! Local variables
704 logical :: are_scalars
705 are_scalars = .false. ; if (present(scalars)) are_scalars = scalars
706
707 if (debug_chksums) then
708 call bchksum(u_comp, mesg//"(u)", g%HI, halos, unscale=unscale)
709 call bchksum(v_comp, mesg//"(v)", g%HI, halos, unscale=unscale)
710 endif
711 if (debug_redundant) then
712 if (are_scalars) then
713 call check_redundant_b(mesg, u_comp, v_comp, g, direction=to_all+scalar_pair, unscale=unscale)
714 else
715 call check_redundant_b(mesg, u_comp, v_comp, g, unscale=unscale)
716 endif
717 endif
718
719end subroutine chksum_vec_b3d
720
721! Do a checksum and redundant point check on a 2d B-grid vector.
722subroutine chksum_vec_b2d(mesg, u_comp, v_comp, G, halos, scalars, symmetric, unscale)
723 character(len=*), intent(in) :: mesg !< An identifying message
724 type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
725 real, dimension(G%IsdB:,G%JsdB:), intent(in) :: u_comp !< The u-component of the vector to be
726 !! checked for consistency in arbitrary,
727 !! possibly rescaled units [A ~> a]
728 real, dimension(G%IsdB:,G%JsdB:), intent(in) :: v_comp !< The v-component of the vector to be
729 !! checked for consistency in arbitrary,
730 !! possibly rescaled units [A ~> a]
731 integer, optional, intent(in) :: halos !< The width of halos to check (default 0)
732 logical, optional, intent(in) :: scalars !< If true this is a pair of
733 !! scalars that are being checked.
734 logical, optional, intent(in) :: symmetric !< If true, do the checksums on the
735 !! full symmetric computational domain.
736 real, optional, intent(in) :: unscale !< A factor that undoes the scaling for the
737 !! arrays to give consistent output [a A-1 ~> 1]
738 ! Local variables
739 logical :: are_scalars
740 are_scalars = .false. ; if (present(scalars)) are_scalars = scalars
741
742 if (debug_chksums) then
743 call bchksum(u_comp, mesg//"(u)", g%HI, halos, symmetric=symmetric, unscale=unscale)
744 call bchksum(v_comp, mesg//"(v)", g%HI, halos, symmetric=symmetric, unscale=unscale)
745 endif
746 if (debug_redundant) then
747 if (are_scalars) then
748 call check_redundant_b(mesg, u_comp, v_comp, g, direction=to_all+scalar_pair, unscale=unscale)
749 else
750 call check_redundant_b(mesg, u_comp, v_comp, g, unscale=unscale)
751 endif
752 endif
753
754end subroutine chksum_vec_b2d
755
756!> Do a checksum and redundant point check on a 3d C-grid vector.
757subroutine chksum_vec_a3d(mesg, u_comp, v_comp, G, halos, scalars, unscale)
758 character(len=*), intent(in) :: mesg !< An identifying message
759 type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
760 real, dimension(G%isd:,G%jsd:,:), intent(in) :: u_comp !< The u-component of the vector to be
761 !! checked for consistency in arbitrary,
762 !! possibly rescaled units [A ~> a]
763 real, dimension(G%isd:,G%jsd:,:), intent(in) :: v_comp !< The v-component of the vector to be
764 !! checked for consistency in arbitrary,
765 !! possibly rescaled units [A ~> a]
766 integer, optional, intent(in) :: halos !< The width of halos to check (default 0)
767 logical, optional, intent(in) :: scalars !< If true this is a pair of
768 !! scalars that are being checked.
769 real, optional, intent(in) :: unscale !< A factor that undoes the scaling for the
770 !! arrays to give consistent output [a A-1 ~> 1]
771 ! Local variables
772 logical :: are_scalars
773 are_scalars = .false. ; if (present(scalars)) are_scalars = scalars
774
775 if (debug_chksums) then
776 call hchksum(u_comp, mesg//"(u)", g%HI, halos, unscale=unscale)
777 call hchksum(v_comp, mesg//"(v)", g%HI, halos, unscale=unscale)
778 endif
779 if (debug_redundant) then
780 if (are_scalars) then
781 call check_redundant_t(mesg, u_comp, v_comp, g, direction=to_all+scalar_pair, unscale=unscale)
782 else
783 call check_redundant_t(mesg, u_comp, v_comp, g, unscale=unscale)
784 endif
785 endif
786
787end subroutine chksum_vec_a3d
788
789!> Do a checksum and redundant point check on a 2d C-grid vector.
790subroutine chksum_vec_a2d(mesg, u_comp, v_comp, G, halos, scalars, unscale)
791 character(len=*), intent(in) :: mesg !< An identifying message
792 type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
793 real, dimension(G%isd:,G%jsd:), intent(in) :: u_comp !< The u-component of the vector to be
794 !! checked for consistency in arbitrary,
795 !! possibly rescaled units [A ~> a]
796 real, dimension(G%isd:,G%jsd:), intent(in) :: v_comp !< The v-component of the vector to be
797 !! checked for consistency in arbitrary,
798 !! possibly rescaled units [A ~> a]
799 integer, optional, intent(in) :: halos !< The width of halos to check (default 0)
800 logical, optional, intent(in) :: scalars !< If true this is a pair of
801 !! scalars that are being checked.
802 real, optional, intent(in) :: unscale !< A factor that undoes the scaling for the
803 !! arrays to give consistent output [a A-1 ~> 1]
804 ! Local variables
805 logical :: are_scalars
806 are_scalars = .false. ; if (present(scalars)) are_scalars = scalars
807
808 if (debug_chksums) then
809 call hchksum(u_comp, mesg//"(u)", g%HI, halos, unscale=unscale)
810 call hchksum(v_comp, mesg//"(v)", g%HI, halos, unscale=unscale)
811 endif
812 if (debug_redundant) then
813 if (are_scalars) then
814 call check_redundant_t(mesg, u_comp, v_comp, g, direction=to_all+scalar_pair, unscale=unscale)
815 else
816 call check_redundant_t(mesg, u_comp, v_comp, g, unscale=unscale)
817 endif
818 endif
819
820end subroutine chksum_vec_a2d
821
822!> This function returns the sum over computational domain of all
823!! processors of hThick*stuff, where stuff is a 3-d array at tracer points.
824function totalstuff(HI, hThick, areaT, stuff, unscale)
825 type(hor_index_type), intent(in) :: hi !< A horizontal index type
826 real, dimension(HI%isd:,HI%jsd:,:), intent(in) :: hthick !< The array of thicknesses to use as weights
827 !! [H ~> m or kg m-2] or [m] or [kg m-2]
828 real, dimension(HI%isd:,HI%jsd:), intent(in) :: areat !< The array of cell areas [L2 ~> m2] or [m2]
829 real, dimension(HI%isd:,HI%jsd:,:), intent(in) :: stuff !< The array of stuff to be summed in arbitrary
830 !! units [A ~> a] or [a]
831 real, optional, intent(in) :: unscale !< A factor that is used to undo scaling of the array
832 !! and the cell mass or volume before it is summed in
833 !! [a m3 A-1 H-1 L-2 ~> 1] or [a kg A-1 H-1 L-2 ~> 1]
834 real :: totalstuff !< the globally integrated amount of stuff
835 !! [A H L2 ~> a m3 or a kg] or [a m3]
836 ! Local variables
837 real :: tmp_for_sum(hi%isc:hi%iec, hi%jsc:hi%jec) ! The column integrated amount of stuff in a
838 ! cell [A H L2 ~> a m3 or a kg] or [a m3]
839 integer :: i, j, k, nz
840
841 nz = size(hthick,3)
842 tmp_for_sum(:,:) = 0.0
843 do k=1,nz ; do j=hi%jsc,hi%jec ; do i=hi%isc,hi%iec
844 tmp_for_sum(i,j) = tmp_for_sum(i,j) + hthick(i,j,k) * stuff(i,j,k) * areat(i,j)
845 enddo ; enddo ; enddo
846 totalstuff = reproducing_sum(tmp_for_sum, unscale=unscale)
847
848end function totalstuff
849
850!> This subroutine display the total thickness, temperature and salinity
851!! as well as the change since the last call.
852subroutine totaltands(HI, hThick, areaT, temperature, salinity, mesg, US, H_to_mks)
853 type(hor_index_type), intent(in) :: hi !< A horizontal index type
854 real, dimension(HI%isd:,HI%jsd:,:), intent(in) :: hthick !< The array of thicknesses to use as weights
855 !! [H ~> m or kg m-2] or [m] or [kg m-2]
856 real, dimension(HI%isd:,HI%jsd:), intent(in) :: areat !< The array of cell areas [L2 ~> m2] or [m2]
857 real, dimension(HI%isd:,HI%jsd:,:), intent(in) :: temperature !< The temperature field to sum [C ~> degC] or [degC]
858 real, dimension(HI%isd:,HI%jsd:,:), intent(in) :: salinity !< The salinity field to sum [S ~> ppt] or [ppt]
859 character(len=*), intent(in) :: mesg !< An identifying message
860 type(unit_scale_type), optional, intent(in) :: us !< A dimensional unit scaling type
861 real, optional, intent(in) :: h_to_mks !< A constant that translates thickness units to its
862 !! MKS units (m or kg m-2) based on whether the model is
863 !! Boussinesq [m H-1 ~> 1] or not [kg m-2 H-1 ~> 1]
864 ! NOTE: This subroutine uses "save" data which is not thread safe and is purely for
865 ! extreme debugging without a proper debugger.
866 real, save :: totalh = 0. ! The total ocean volume or mass, saved for the next
867 ! call [H L2 ~> m3 or kg] or [m3] or [kg]
868 real, save :: totalt = 0. ! The total volume integrated ocean temperature, saved for the next
869 ! call [C H L2 ~> degC m3 or degC kg] or [degC m3] or [degC kg]
870 real, save :: totals = 0. ! The total volume integrated ocean salinity, saved for the next
871 ! call [S H L2 ~> ppt m3 or ppt kg] or [ppt m3] or [ppt kg]
872 ! Local variables
873 logical, save :: firstcall = .true.
874 real :: tmp_for_sum(hi%isc:hi%iec, hi%jsc:hi%jec) ! The volume of each column [H L2 ~> m3 or kg] or [m3] or [kg]
875 real :: thish, delh ! The total ocean volume and the change from the last call [H L2 ~> m3 or kg] or [m3] or [kg]
876 real :: thist, delt ! The current total volume integrated temperature and the change from the last
877 ! call [C H L2 ~> degC m3 or degC kg] or [degC m3] or [degC kg]
878 real :: thiss, dels ! The current total volume integrated salinity and the change from the last
879 ! call [S H L2 ~> ppt m3 or ppt kg] or [ppt m3] or [ppt kg]
880 real :: h_unscale ! A constant that translates thickness units to its MKS units (m or kg m-2) based on
881 ! whether the model is Boussinesq [m H-1 ~> 1] or non-Boussinesq [kg m-2 H-1 ~> 1]
882 real :: hl2_unscale ! An overall unscaling factor for cell mass or volume [m3 H-1 L-2 ~> 1] or [kg H-1 L-2 ~> 1]
883 real :: t_unscale ! An overall unscaling factor for cell-integrated temperature [degC m3 C-1 H-1 L-2 ~> 1] or
884 ! [degC kg C-1 H-1 L-2 ~> 1]
885 real :: s_unscale ! An overall unscaling factor for cell-integrated salinity [ppt m3 S-1 H-1 L-2 ~> 1] or
886 ! [ppt kg S-1 H-1 L-2 ~> 1]
887 integer :: i, j, k, nz
888
889 h_unscale = 1.0 ; if (present(h_to_mks)) h_unscale = h_to_mks
890 if (present(us)) then
891 hl2_unscale = us%L_to_m**2 * h_unscale
892 t_unscale = us%C_to_degC * hl2_unscale ; s_unscale = us%S_to_ppt * hl2_unscale
893 else
894 hl2_unscale = h_unscale
895 t_unscale = hl2_unscale ; s_unscale = hl2_unscale
896 endif
897
898 nz = size(hthick,3)
899 tmp_for_sum(:,:) = 0.0
900 do k=1,nz ; do j=hi%jsc,hi%jec ; do i=hi%isc,hi%iec
901 tmp_for_sum(i,j) = tmp_for_sum(i,j) + hthick(i,j,k) * areat(i,j)
902 enddo ; enddo ; enddo
903 thish = reproducing_sum(tmp_for_sum, unscale=hl2_unscale)
904 thist = totalstuff(hi, hthick, areat, temperature, unscale=t_unscale)
905 thiss = totalstuff(hi, hthick, areat, salinity, unscale=s_unscale)
906
907 if (is_root_pe()) then
908 if (firstcall) then
909 totalh = thish ; totalt = thist ; totals = thiss
910 write(stdout,*) 'Totals H,T,S:', thish*hl2_unscale, thist*t_unscale, thiss*s_unscale, ' ', mesg
911 firstcall = .false.
912 else
913 delh = thish - totalh
914 delt = thist - totalt
915 dels = thiss - totals
916 totalh = thish ; totalt = thist ; totals = thiss
917 write(0,*) 'Tot/del H,T,S:', thish*hl2_unscale, thist*t_unscale, thiss*s_unscale, &
918 delh*hl2_unscale, delt*t_unscale, dels*s_unscale, ' ', mesg
919 endif
920 endif
921
922end subroutine totaltands
923
924!> Returns false if the column integral of a given quantity is within roundoff
925logical function check_column_integral(nk, field, known_answer)
926 integer, intent(in) :: nk !< Number of levels in column
927 real, dimension(nk), intent(in) :: field !< Field to be summed [arbitrary]
928 real, optional, intent(in) :: known_answer !< If present is the expected sum [arbitrary],
929 !! If missing, assumed zero
930 ! Local variables
931 real :: u_sum ! The vertical sum of the field [arbitrary]
932 real :: error ! An estimate of the roundoff error in the sum [arbitrary]
933 real :: expected ! The expected vertical sum [arbitrary]
934 integer :: k
935
936 u_sum = field(1)
937 error = 0.
938
939 ! Reintegrate and sum roundoff errors
940 do k=2,nk
941 u_sum = u_sum + field(k)
942 error = error + epsilon(u_sum)*max(abs(u_sum),abs(field(k)))
943 enddo
944
945 ! Assign expected answer to either the optional input or 0
946 if (present(known_answer)) then
947 expected = known_answer
948 else
949 expected = 0.
950 endif
951
952 ! Compare the column integrals against calculated roundoff error
953 if (abs(u_sum-expected) > error) then
954 check_column_integral = .true.
955 else
956 check_column_integral = .false.
957 endif
958
959end function check_column_integral
960
961!> Returns false if the column integrals of two given quantities are within roundoff of each other
962logical function check_column_integrals(nk_1, field_1, nk_2, field_2, missing_value)
963 integer, intent(in) :: nk_1 !< Number of levels in field 1
964 integer, intent(in) :: nk_2 !< Number of levels in field 2
965 real, dimension(nk_1), intent(in) :: field_1 !< First field to be summed [arbitrary]
966 real, dimension(nk_2), intent(in) :: field_2 !< Second field to be summed [arbitrary]
967 real, optional, intent(in) :: missing_value !< If column contains missing values,
968 !! mask them from the sum [arbitrary]
969 ! Local variables
970 real :: u1_sum, u2_sum ! The vertical sums of the two fields [arbitrary]
971 real :: error1, error2 ! Estimates of the roundoff errors in the sums [arbitrary]
972 real :: misval ! The missing value flag, indicating elements that are to be omitted
973 ! from the sums [arbitrary]
974 integer :: k
975
976 ! Assign missing value
977 if (present(missing_value)) then
978 misval = missing_value
979 else
980 misval = 0.
981 endif
982
983 u1_sum = field_1(1)
984 error1 = 0.
985
986 ! Reintegrate and sum roundoff errors
987 do k=2,nk_1
988 if (field_1(k) /= misval) then
989 u1_sum = u1_sum + field_1(k)
990 error1 = error1 + epsilon(u1_sum)*max(abs(u1_sum),abs(field_1(k)))
991 endif
992 enddo
993
994 u2_sum = field_2(1)
995 error2 = 0.
996
997 ! Reintegrate and sum roundoff errors
998 do k=2,nk_2
999 if (field_2(k) /= misval) then
1000 u2_sum = u2_sum + field_2(k)
1001 error2 = error2 + epsilon(u2_sum)*max(abs(u2_sum),abs(field_2(k)))
1002 endif
1003 enddo
1004
1005 ! Compare the column integrals against calculated roundoff error
1006 if (abs(u1_sum-u2_sum) > (error1+error2)) then
1007 check_column_integrals = .true.
1008 else
1009 check_column_integrals = .false.
1010 endif
1011
1012end function check_column_integrals
1013
1014end module mom_debugging