portedportable, not yet portedexecuted, not portableexecutable, not hit by this run
| 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. | |
| 11 | module MOM_debugging | |
| 12 | ||
| 13 | use MOM_checksums, only : hchksum, Bchksum, qchksum, uvchksum, hchksum_pair | |
| 14 | use MOM_checksums, only : is_NaN, chksum, MOM_checksums_init | |
| 15 | use MOM_coms, only : PE_here, root_PE, num_PEs | |
| 16 | use MOM_coms, only : min_across_PEs, max_across_PEs, reproducing_sum | |
| 17 | use MOM_domains, only : pass_vector, pass_var, pe_here | |
| 18 | use MOM_domains, only : BGRID_NE, AGRID, To_All, Scalar_Pair | |
| 19 | use MOM_error_handler, only : MOM_error, FATAL, WARNING, is_root_pe | |
| 20 | use MOM_file_parser, only : log_version, param_file_type, get_param | |
| 21 | use MOM_grid, only : ocean_grid_type | |
| 22 | use MOM_hor_index, only : hor_index_type | |
| 23 | use MOM_io, only : stdout | |
| 24 | use MOM_unit_scaling, only : unit_scale_type | |
| 25 | ||
| 26 | implicit none ; private | |
| 27 | ||
| 28 | public :: check_redundant_C, check_redundant_B, check_redundant_T, check_redundant | |
| 29 | public :: vec_chksum, vec_chksum_C, vec_chksum_B, vec_chksum_A | |
| 30 | public :: MOM_debugging_init, totalStuff, totalTandS | |
| 31 | public :: check_column_integral, check_column_integrals | |
| 32 | public :: query_debugging_checks | |
| 33 | ||
| 34 | ! These interfaces come from MOM_checksums. | |
| 35 | public :: hchksum, Bchksum, qchksum, is_NaN, chksum, uvchksum, hchksum_pair | |
| 36 | ||
| 37 | !> Check for consistency between the duplicated points of a C-grid vector | |
| 38 | interface check_redundant | |
| 39 | module procedure check_redundant_vC3d, check_redundant_vC2d | |
| 40 | end interface check_redundant | |
| 41 | !> Check for consistency between the duplicated points of a C-grid vector | |
| 42 | interface check_redundant_C | |
| 43 | module procedure check_redundant_vC3d, check_redundant_vC2d | |
| 44 | end interface check_redundant_C | |
| 45 | !> Check for consistency between the duplicated points of a B-grid vector or scalar | |
| 46 | interface check_redundant_B | |
| 47 | module procedure check_redundant_vB3d, check_redundant_vB2d | |
| 48 | module procedure check_redundant_sB3d, check_redundant_sB2d | |
| 49 | end interface check_redundant_B | |
| 50 | !> Check for consistency between the duplicated points of an A-grid vector or scalar | |
| 51 | interface check_redundant_T | |
| 52 | module procedure check_redundant_sT3d, check_redundant_sT2d | |
| 53 | module procedure check_redundant_vT3d, check_redundant_vT2d | |
| 54 | end interface check_redundant_T | |
| 55 | ||
| 56 | !> Do checksums on the components of a C-grid vector | |
| 57 | interface vec_chksum | |
| 58 | module procedure chksum_vec_C3d, chksum_vec_C2d | |
| 59 | end interface vec_chksum | |
| 60 | !> Do checksums on the components of a C-grid vector | |
| 61 | interface vec_chksum_C | |
| 62 | module procedure chksum_vec_C3d, chksum_vec_C2d | |
| 63 | end interface vec_chksum_C | |
| 64 | !> Do checksums on the components of a B-grid vector | |
| 65 | interface vec_chksum_B | |
| 66 | module procedure chksum_vec_B3d, chksum_vec_B2d | |
| 67 | end interface vec_chksum_B | |
| 68 | !> Do checksums on the components of an A-grid vector | |
| 69 | interface vec_chksum_A | |
| 70 | module procedure chksum_vec_A3d, chksum_vec_A2d | |
| 71 | end 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. | |
| 75 | integer :: max_redundant_prints = 100 !< Maximum number of times to write redundant messages | |
| 76 | integer :: redundant_prints(3) = 0 !< Counters for controlling redundant printing | |
| 77 | logical :: debug = .false. !< Write out verbose debugging data | |
| 78 | logical :: debug_chksums = .true. !< Perform checksums on arrays | |
| 79 | logical :: debug_redundant = .true. !< Check redundant values on PE boundaries | |
| 80 | ||
| 81 | contains | |
| 82 | ||
| 83 | !> MOM_debugging_init initializes the MOM_debugging module, and sets | |
| 84 | !! the parameters that control which checks are active for MOM6. | |
| 85 | 1 | subroutine 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 | 1 | 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 | 1 | 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 | 1 | 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 | 1 | debuggingParam=.true.) |
| 103 | ||
| 104 | 1 | call MOM_checksums_init(param_file) |
| 105 | ||
| 106 | 1 | end subroutine MOM_debugging_init |
| 107 | ||
| 108 | !> Returns logicals indicating which debugging checks should be performed. | |
| 109 | 0 | subroutine 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 | 0 | if (present(do_debug)) do_debug = debug |
| 115 | 0 | if (present(do_chksums)) do_chksums = debug_chksums |
| 116 | 0 | if (present(do_redundant)) do_redundant = debug_redundant |
| 117 | ||
| 118 | 0 | end subroutine query_debugging_checks |
| 119 | ||
| 120 | !> Check for consistency between the duplicated points of a 3-D C-grid vector | |
| 121 | 0 | subroutine 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 | 0 | do k=1,size(u_comp,3) |
| 145 | 0 | write(mesg_k,'(" Layer ",i0," ")') k |
| 146 | call check_redundant_vC2d(trim(mesg)//trim(mesg_k), u_comp(:,:,k), & | |
| 147 | 0 | v_comp(:,:,k), G, is, ie, js, je, direction, unscale) |
| 148 | enddo | |
| 149 | 0 | end subroutine check_redundant_vC3d |
| 150 | ||
| 151 | !> Check for consistency between the duplicated points of a 2-D C-grid vector | |
| 152 | 0 | subroutine 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 | 0 | real :: u_nonsym(G%isd:G%ied,G%jsd:G%jed) ! A nonsymmetric version of u_comp [A ~> a] |
| 174 | 0 | real :: v_nonsym(G%isd:G%ied,G%jsd:G%jed) ! A nonsymmetric version of v_comp [A ~> a] |
| 175 | 0 | real :: u_resym(G%IsdB:G%IedB,G%jsd:G%jed) ! A reconstructed symmetric version of u_comp [A ~> a] |
| 176 | 0 | 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 | 0 | Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB |
| 183 | 0 | isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed |
| 184 | 0 | IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB |
| 185 | ||
| 186 | 0 | 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 | 0 | if ((isd == IsdB) .and. (jsd == JsdB)) return |
| 189 | endif | |
| 190 | ||
| 191 | 0 | sc = 1.0 ; if (present(unscale)) sc = unscale |
| 192 | ||
| 193 | 0 | do i=isd,ied ; do j=jsd,jed |
| 194 | 0 | u_nonsym(i,j) = u_comp(i,j) ; v_nonsym(i,j) = v_comp(i,j) |
| 195 | enddo ; enddo | |
| 196 | ||
| 197 | 0 | if (.not.associated(G%Domain_aux)) call MOM_error(FATAL, & |
| 198 | 0 | " check_redundant called with a non-associated auxiliary domain the grid type.") |
| 199 | 0 | call pass_vector(u_nonsym, v_nonsym, G%Domain_aux, direction) |
| 200 | ||
| 201 | 0 | do I=IsdB,IedB ; do j=jsd,jed ; u_resym(I,j) = u_comp(I,j) ; enddo ; enddo |
| 202 | 0 | do i=isd,ied ; do J=JsdB,JedB ; v_resym(i,J) = v_comp(i,J) ; enddo ; enddo |
| 203 | 0 | do i=isd,ied ; do j=jsd,jed |
| 204 | 0 | u_resym(i,j) = u_nonsym(i,j) ; v_resym(i,j) = v_nonsym(i,j) |
| 205 | enddo ; enddo | |
| 206 | 0 | call pass_vector(u_resym, v_resym, G%Domain, direction) |
| 207 | ||
| 208 | 0 | is_ch = Isq ; ie_ch = Ieq ; js_ch = Jsq ; je_ch = Jeq |
| 209 | 0 | if (present(is)) is_ch = is ; if (present(ie)) ie_ch = ie |
| 210 | 0 | if (present(js)) js_ch = js ; if (present(js)) je_ch = je |
| 211 | ||
| 212 | 0 | do i=is_ch,ie_ch ; do j=js_ch+1,je_ch |
| 213 | 0 | if (u_resym(i,j) /= u_comp(i,j) .and. & |
| 214 | 0 | 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 | 0 | sc*u_comp(i,j), sc*u_resym(i,j), sc*(u_comp(i,j)-u_resym(i,j)), i, j, pe_here() |
| 218 | 0 | write(0,'(A130)') trim(mesg)//trim(mesg2) |
| 219 | 0 | redundant_prints(3) = redundant_prints(3) + 1 |
| 220 | endif | |
| 221 | enddo ; enddo | |
| 222 | 0 | do i=is_ch+1,ie_ch ; do j=js_ch,je_ch |
| 223 | 0 | if (v_resym(i,j) /= v_comp(i,j) .and. & |
| 224 | 0 | 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 | 0 | sc*v_comp(i,j), sc*v_resym(i,j), sc*(v_comp(i,j)-v_resym(i,j)), i, j, & |
| 228 | 0 | G%geoLonBu(i,j), G%geoLatBu(i,j), pe_here() |
| 229 | 0 | write(0,'(A155)') trim(mesg)//trim(mesg2) |
| 230 | 0 | redundant_prints(3) = redundant_prints(3) + 1 |
| 231 | endif | |
| 232 | enddo ; enddo | |
| 233 | ||
| 234 | 0 | end subroutine check_redundant_vC2d |
| 235 | ||
| 236 | !> Check for consistency between the duplicated points of a 3-D scalar at corner points | |
| 237 | 0 | subroutine 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 | 0 | do k=1,size(array,3) |
| 254 | 0 | write(mesg_k,'(" Layer ",i0," ")') k |
| 255 | call check_redundant_sB2d(trim(mesg)//trim(mesg_k), array(:,:,k), & | |
| 256 | 0 | G, is, ie, js, je, unscale) |
| 257 | enddo | |
| 258 | 0 | end subroutine check_redundant_sB3d |
| 259 | ||
| 260 | !> Check for consistency between the duplicated points of a 2-D scalar at corner points | |
| 261 | 0 | subroutine 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 | 0 | real :: a_nonsym(G%isd:G%ied,G%jsd:G%jed) ! A nonsymmetric version of array [A ~> a] |
| 276 | 0 | 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 | 0 | Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB |
| 283 | 0 | isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed |
| 284 | 0 | IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB |
| 285 | ||
| 286 | 0 | 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 | 0 | if ((isd == IsdB) .and. (jsd == JsdB)) return |
| 289 | endif | |
| 290 | ||
| 291 | 0 | sc = 1.0 ; if (present(unscale)) sc = unscale |
| 292 | ||
| 293 | 0 | do i=isd,ied ; do j=jsd,jed |
| 294 | 0 | a_nonsym(i,j) = array(i,j) |
| 295 | enddo ; enddo | |
| 296 | ||
| 297 | 0 | if (.not.associated(G%Domain_aux)) call MOM_error(FATAL, & |
| 298 | 0 | " check_redundant called with a non-associated auxiliary domain the grid type.") |
| 299 | call pass_vector(a_nonsym, a_nonsym, G%Domain_aux, & | |
| 300 | 0 | direction=To_All+Scalar_Pair, stagger=BGRID_NE) |
| 301 | ||
| 302 | 0 | do I=IsdB,IedB ; do J=JsdB,JedB ; a_resym(I,J) = array(I,J) ; enddo ; enddo |
| 303 | 0 | do i=isd,ied ; do j=jsd,jed |
| 304 | 0 | 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 | 0 | stagger=BGRID_NE) |
| 308 | ||
| 309 | 0 | is_ch = Isq ; ie_ch = Ieq ; js_ch = Jsq ; je_ch = Jeq |
| 310 | 0 | if (present(is)) is_ch = is ; if (present(ie)) ie_ch = ie |
| 311 | 0 | if (present(js)) js_ch = js ; if (present(js)) je_ch = je |
| 312 | ||
| 313 | 0 | do i=is_ch,ie_ch ; do j=js_ch,je_ch |
| 314 | 0 | if (a_resym(i,j) /= array(i,j) .and. & |
| 315 | 0 | 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 | 0 | sc*array(i,j), sc*a_resym(i,j), sc*(array(i,j)-a_resym(i,j)), i, j, pe_here() |
| 319 | 0 | write(0,'(A130)') trim(mesg)//trim(mesg2) |
| 320 | 0 | redundant_prints(2) = redundant_prints(2) + 1 |
| 321 | endif | |
| 322 | enddo ; enddo | |
| 323 | ||
| 324 | 0 | end subroutine check_redundant_sB2d |
| 325 | ||
| 326 | !> Check for consistency between the duplicated points of a 3-D B-grid vector | |
| 327 | 0 | subroutine 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 | 0 | do k=1,size(u_comp,3) |
| 350 | 0 | write(mesg_k,'(" Layer ",i0," ")') k |
| 351 | call check_redundant_vB2d(trim(mesg)//trim(mesg_k), u_comp(:,:,k), & | |
| 352 | 0 | v_comp(:,:,k), G, is, ie, js, je, direction, unscale) |
| 353 | enddo | |
| 354 | 0 | end subroutine check_redundant_vB3d |
| 355 | ||
| 356 | !> Check for consistency between the duplicated points of a 2-D B-grid vector | |
| 357 | 0 | subroutine 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 | 0 | real :: u_nonsym(G%isd:G%ied,G%jsd:G%jed) ! A nonsymmetric version of u_comp [A ~> a] |
| 379 | 0 | real :: v_nonsym(G%isd:G%ied,G%jsd:G%jed) ! A nonsymmetric version of v_comp [A ~> a] |
| 380 | 0 | real :: u_resym(G%IsdB:G%IedB,G%JsdB:G%JedB) ! A reconstructed symmetric version of u_comp [A ~> a] |
| 381 | 0 | 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 | 0 | Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB |
| 388 | 0 | isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed |
| 389 | 0 | IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB |
| 390 | ||
| 391 | 0 | 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 | 0 | if ((isd == IsdB) .and. (jsd == JsdB)) return |
| 394 | endif | |
| 395 | ||
| 396 | 0 | sc = 1.0 ; if (present(unscale)) sc = unscale |
| 397 | ||
| 398 | 0 | do i=isd,ied ; do j=jsd,jed |
| 399 | 0 | u_nonsym(i,j) = u_comp(i,j) ; v_nonsym(i,j) = v_comp(i,j) |
| 400 | enddo ; enddo | |
| 401 | ||
| 402 | 0 | if (.not.associated(G%Domain_aux)) call MOM_error(FATAL, & |
| 403 | 0 | " check_redundant called with a non-associated auxiliary domain the grid type.") |
| 404 | 0 | call pass_vector(u_nonsym, v_nonsym, G%Domain_aux, direction, stagger=BGRID_NE) |
| 405 | ||
| 406 | 0 | do I=IsdB,IedB ; do J=JsdB,JedB |
| 407 | 0 | u_resym(I,J) = u_comp(I,J) ; v_resym(I,J) = v_comp(I,J) |
| 408 | enddo ; enddo | |
| 409 | 0 | do i=isd,ied ; do j=jsd,jed |
| 410 | 0 | u_resym(i,j) = u_nonsym(i,j) ; v_resym(i,j) = v_nonsym(i,j) |
| 411 | enddo ; enddo | |
| 412 | 0 | call pass_vector(u_resym, v_resym, G%Domain, direction, stagger=BGRID_NE) |
| 413 | ||
| 414 | 0 | is_ch = Isq ; ie_ch = Ieq ; js_ch = Jsq ; je_ch = Jeq |
| 415 | 0 | if (present(is)) is_ch = is ; if (present(ie)) ie_ch = ie |
| 416 | 0 | if (present(js)) js_ch = js ; if (present(js)) je_ch = je |
| 417 | ||
| 418 | 0 | do i=is_ch,ie_ch ; do j=js_ch,je_ch |
| 419 | 0 | if (u_resym(i,j) /= u_comp(i,j) .and. & |
| 420 | 0 | 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 | 0 | sc*u_comp(i,j), sc*u_resym(i,j), sc*(u_comp(i,j)-u_resym(i,j)), i, j, pe_here() |
| 424 | 0 | write(0,'(A130)') trim(mesg)//trim(mesg2) |
| 425 | 0 | redundant_prints(2) = redundant_prints(2) + 1 |
| 426 | endif | |
| 427 | enddo ; enddo | |
| 428 | 0 | do i=is_ch,ie_ch ; do j=js_ch,je_ch |
| 429 | 0 | if (v_resym(i,j) /= v_comp(i,j) .and. & |
| 430 | 0 | 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 | 0 | sc*v_comp(i,j), sc*v_resym(i,j), sc*(v_comp(i,j)-v_resym(i,j)), i, j, & |
| 434 | 0 | G%geoLonBu(i,j), G%geoLatBu(i,j), pe_here() |
| 435 | 0 | write(0,'(A155)') trim(mesg)//trim(mesg2) |
| 436 | 0 | redundant_prints(2) = redundant_prints(2) + 1 |
| 437 | endif | |
| 438 | enddo ; enddo | |
| 439 | ||
| 440 | 0 | end subroutine check_redundant_vB2d |
| 441 | ||
| 442 | !> Check for consistency between the duplicated points of a 3-D scalar at tracer points | |
| 443 | 0 | subroutine 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 | 0 | do k=1,size(array,3) |
| 459 | 0 | write(mesg_k,'(" Layer ",i0," ")') k |
| 460 | call check_redundant_sT2d(trim(mesg)//trim(mesg_k), array(:,:,k), & | |
| 461 | 0 | G, is, ie, js, je, unscale) |
| 462 | enddo | |
| 463 | 0 | end subroutine check_redundant_sT3d |
| 464 | ||
| 465 | ||
| 466 | !> Check for consistency between the duplicated points of a 2-D scalar at tracer points | |
| 467 | 0 | subroutine 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 | 0 | 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 | 0 | isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed |
| 488 | ||
| 489 | 0 | is_ch = G%isc ; ie_ch = G%iec ; js_ch = G%jsc ; je_ch = G%jec |
| 490 | 0 | if (present(is)) is_ch = is ; if (present(ie)) ie_ch = ie |
| 491 | 0 | if (present(js)) js_ch = js ; if (present(js)) je_ch = je |
| 492 | ||
| 493 | 0 | 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 | 0 | (js_ch == G%jsc) .and. (je_ch == G%jec)) return |
| 498 | ||
| 499 | 0 | do i=isd,ied ; do j=jsd,jed |
| 500 | 0 | a_nonsym(i,j) = array(i,j) |
| 501 | enddo ; enddo | |
| 502 | ||
| 503 | 0 | call pass_var(a_nonsym, G%Domain) |
| 504 | ||
| 505 | 0 | do i=is_ch,ie_ch ; do j=js_ch,je_ch |
| 506 | 0 | if (a_nonsym(i,j) /= array(i,j) .and. & |
| 507 | 0 | 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 | 0 | sc*array(i,j), sc*a_nonsym(i,j), sc*(array(i,j)-a_nonsym(i,j)), i, j, pe_here() |
| 511 | 0 | write(0,'(A130)') trim(mesg)//trim(mesg2) |
| 512 | 0 | redundant_prints(1) = redundant_prints(1) + 1 |
| 513 | endif | |
| 514 | enddo ; enddo | |
| 515 | ||
| 516 | 0 | end subroutine check_redundant_sT2d |
| 517 | ||
| 518 | !> Check for consistency between the duplicated points of a 3-D A-grid vector | |
| 519 | 0 | subroutine 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 | 0 | do k=1,size(u_comp,3) |
| 542 | 0 | write(mesg_k,'(" Layer ",i0," ")') k |
| 543 | call check_redundant_vT2d(trim(mesg)//trim(mesg_k), u_comp(:,:,k), & | |
| 544 | 0 | v_comp(:,:,k), G, is, ie, js, je, direction, unscale) |
| 545 | enddo | |
| 546 | 0 | end subroutine check_redundant_vT3d |
| 547 | ||
| 548 | !> Check for consistency between the duplicated points of a 2-D A-grid vector | |
| 549 | 0 | subroutine 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 | 0 | 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 | 0 | 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 | 0 | Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB |
| 578 | 0 | isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed |
| 579 | 0 | IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB |
| 580 | ||
| 581 | 0 | is_ch = G%isc ; ie_ch = G%iec ; js_ch = G%jsc ; je_ch = G%jec |
| 582 | 0 | if (present(is)) is_ch = is ; if (present(ie)) ie_ch = ie |
| 583 | 0 | if (present(js)) js_ch = js ; if (present(js)) je_ch = je |
| 584 | ||
| 585 | 0 | 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 | 0 | (js_ch == G%jsc) .and. (je_ch == G%jec)) return |
| 590 | ||
| 591 | 0 | do i=isd,ied ; do j=jsd,jed |
| 592 | 0 | u_nonsym(i,j) = u_comp(i,j) ; v_nonsym(i,j) = v_comp(i,j) |
| 593 | enddo ; enddo | |
| 594 | ||
| 595 | 0 | call pass_vector(u_nonsym, v_nonsym, G%Domain, direction, stagger=AGRID) |
| 596 | ||
| 597 | 0 | do i=is_ch,ie_ch ; do j=js_ch+1,je_ch |
| 598 | 0 | if (u_nonsym(i,j) /= u_comp(i,j) .and. & |
| 599 | 0 | 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 | 0 | sc*u_comp(i,j), sc*u_nonsym(i,j), sc*(u_comp(i,j)-u_nonsym(i,j)), i, j, pe_here() |
| 603 | 0 | write(0,'(A130)') trim(mesg)//trim(mesg2) |
| 604 | 0 | redundant_prints(1) = redundant_prints(1) + 1 |
| 605 | endif | |
| 606 | enddo ; enddo | |
| 607 | 0 | do i=is_ch+1,ie_ch ; do j=js_ch,je_ch |
| 608 | 0 | if (v_nonsym(i,j) /= v_comp(i,j) .and. & |
| 609 | 0 | 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 | 0 | sc*v_comp(i,j), sc*v_nonsym(i,j), sc*(v_comp(i,j)-v_nonsym(i,j)), i, j, & |
| 613 | 0 | G%geoLonBu(i,j), G%geoLatBu(i,j), pe_here() |
| 614 | 0 | write(0,'(A155)') trim(mesg)//trim(mesg2) |
| 615 | 0 | redundant_prints(1) = redundant_prints(1) + 1 |
| 616 | endif | |
| 617 | enddo ; enddo | |
| 618 | ||
| 619 | 0 | end 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. | |
| 625 | 0 | subroutine 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 | 0 | are_scalars = .false. ; if (present(scalars)) are_scalars = scalars |
| 642 | ||
| 643 | 0 | if (debug_chksums) then |
| 644 | 0 | call uvchksum(mesg, u_comp, v_comp, G%HI, halos, unscale=unscale) |
| 645 | endif | |
| 646 | 0 | if (debug_redundant) then |
| 647 | 0 | if (are_scalars) then |
| 648 | 0 | call check_redundant_C(mesg, u_comp, v_comp, G, direction=To_All+Scalar_Pair, unscale=unscale) |
| 649 | else | |
| 650 | 0 | call check_redundant_C(mesg, u_comp, v_comp, G, unscale=unscale) |
| 651 | endif | |
| 652 | endif | |
| 653 | ||
| 654 | 0 | end subroutine chksum_vec_C3d |
| 655 | ||
| 656 | !> Do a checksum and redundant point check on a 2d C-grid vector. | |
| 657 | 0 | subroutine 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 | 0 | are_scalars = .false. ; if (present(scalars)) are_scalars = scalars |
| 674 | ||
| 675 | 0 | if (debug_chksums) then |
| 676 | 0 | call uvchksum(mesg, u_comp, v_comp, G%HI, halos, unscale=unscale) |
| 677 | endif | |
| 678 | 0 | if (debug_redundant) then |
| 679 | 0 | if (are_scalars) then |
| 680 | 0 | call check_redundant_C(mesg, u_comp, v_comp, G, direction=To_All+Scalar_Pair, unscale=unscale) |
| 681 | else | |
| 682 | 0 | call check_redundant_C(mesg, u_comp, v_comp, G, unscale=unscale) |
| 683 | endif | |
| 684 | endif | |
| 685 | ||
| 686 | 0 | end subroutine chksum_vec_C2d |
| 687 | ||
| 688 | !> Do a checksum and redundant point check on a 3d B-grid vector. | |
| 689 | 0 | subroutine 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 | 0 | are_scalars = .false. ; if (present(scalars)) are_scalars = scalars |
| 706 | ||
| 707 | 0 | if (debug_chksums) then |
| 708 | 0 | call Bchksum(u_comp, mesg//"(u)", G%HI, halos, unscale=unscale) |
| 709 | 0 | call Bchksum(v_comp, mesg//"(v)", G%HI, halos, unscale=unscale) |
| 710 | endif | |
| 711 | 0 | if (debug_redundant) then |
| 712 | 0 | if (are_scalars) then |
| 713 | 0 | call check_redundant_B(mesg, u_comp, v_comp, G, direction=To_All+Scalar_Pair, unscale=unscale) |
| 714 | else | |
| 715 | 0 | call check_redundant_B(mesg, u_comp, v_comp, G, unscale=unscale) |
| 716 | endif | |
| 717 | endif | |
| 718 | ||
| 719 | 0 | end subroutine chksum_vec_B3d |
| 720 | ||
| 721 | ! Do a checksum and redundant point check on a 2d B-grid vector. | |
| 722 | 0 | subroutine 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 | 0 | are_scalars = .false. ; if (present(scalars)) are_scalars = scalars |
| 741 | ||
| 742 | 0 | if (debug_chksums) then |
| 743 | 0 | call Bchksum(u_comp, mesg//"(u)", G%HI, halos, symmetric=symmetric, unscale=unscale) |
| 744 | 0 | call Bchksum(v_comp, mesg//"(v)", G%HI, halos, symmetric=symmetric, unscale=unscale) |
| 745 | endif | |
| 746 | 0 | if (debug_redundant) then |
| 747 | 0 | if (are_scalars) then |
| 748 | 0 | call check_redundant_B(mesg, u_comp, v_comp, G, direction=To_All+Scalar_Pair, unscale=unscale) |
| 749 | else | |
| 750 | 0 | call check_redundant_B(mesg, u_comp, v_comp, G, unscale=unscale) |
| 751 | endif | |
| 752 | endif | |
| 753 | ||
| 754 | 0 | end subroutine chksum_vec_B2d |
| 755 | ||
| 756 | !> Do a checksum and redundant point check on a 3d C-grid vector. | |
| 757 | 0 | subroutine 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 | 0 | are_scalars = .false. ; if (present(scalars)) are_scalars = scalars |
| 774 | ||
| 775 | 0 | if (debug_chksums) then |
| 776 | 0 | call hchksum(u_comp, mesg//"(u)", G%HI, halos, unscale=unscale) |
| 777 | 0 | call hchksum(v_comp, mesg//"(v)", G%HI, halos, unscale=unscale) |
| 778 | endif | |
| 779 | 0 | if (debug_redundant) then |
| 780 | 0 | if (are_scalars) then |
| 781 | 0 | call check_redundant_T(mesg, u_comp, v_comp, G, direction=To_All+Scalar_Pair, unscale=unscale) |
| 782 | else | |
| 783 | 0 | call check_redundant_T(mesg, u_comp, v_comp, G, unscale=unscale) |
| 784 | endif | |
| 785 | endif | |
| 786 | ||
| 787 | 0 | end subroutine chksum_vec_A3d |
| 788 | ||
| 789 | !> Do a checksum and redundant point check on a 2d C-grid vector. | |
| 790 | 0 | subroutine 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 | 0 | are_scalars = .false. ; if (present(scalars)) are_scalars = scalars |
| 807 | ||
| 808 | 0 | if (debug_chksums) then |
| 809 | 0 | call hchksum(u_comp, mesg//"(u)", G%HI, halos, unscale=unscale) |
| 810 | 0 | call hchksum(v_comp, mesg//"(v)", G%HI, halos, unscale=unscale) |
| 811 | endif | |
| 812 | 0 | if (debug_redundant) then |
| 813 | 0 | if (are_scalars) then |
| 814 | 0 | call check_redundant_T(mesg, u_comp, v_comp, G, direction=To_All+Scalar_Pair, unscale=unscale) |
| 815 | else | |
| 816 | 0 | call check_redundant_T(mesg, u_comp, v_comp, G, unscale=unscale) |
| 817 | endif | |
| 818 | endif | |
| 819 | ||
| 820 | 0 | end 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. | |
| 824 | 0 | function 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 | 0 | 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 | 0 | nz = size(hThick,3) |
| 842 | 0 | tmp_for_sum(:,:) = 0.0 |
| 843 | 0 | do k=1,nz ; do j=HI%jsc,HI%jec ; do i=HI%isc,HI%iec |
| 844 | 0 | 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 | 0 | totalStuff = reproducing_sum(tmp_for_sum, unscale=unscale) |
| 847 | ||
| 848 | 0 | end function totalStuff |
| 849 | ||
| 850 | !> This subroutine display the total thickness, temperature and salinity | |
| 851 | !! as well as the change since the last call. | |
| 852 | 0 | subroutine 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 | 0 | 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 | 0 | H_unscale = 1.0 ; if (present(H_to_mks)) H_unscale = H_to_mks |
| 890 | 0 | if (present(US)) then |
| 891 | 0 | HL2_unscale = US%L_to_m**2 * H_unscale |
| 892 | 0 | T_unscale = US%C_to_degC * HL2_unscale ; S_unscale = US%S_to_ppt * HL2_unscale |
| 893 | else | |
| 894 | 0 | HL2_unscale = H_unscale |
| 895 | 0 | T_unscale = HL2_unscale ; S_unscale = HL2_unscale |
| 896 | endif | |
| 897 | ||
| 898 | 0 | nz = size(hThick,3) |
| 899 | 0 | tmp_for_sum(:,:) = 0.0 |
| 900 | 0 | do k=1,nz ; do j=HI%jsc,HI%jec ; do i=HI%isc,HI%iec |
| 901 | 0 | tmp_for_sum(i,j) = tmp_for_sum(i,j) + hThick(i,j,k) * areaT(i,j) |
| 902 | enddo ; enddo ; enddo | |
| 903 | 0 | thisH = reproducing_sum(tmp_for_sum, unscale=HL2_unscale) |
| 904 | 0 | thisT = totalStuff(HI, hThick, areaT, temperature, unscale=T_unscale) |
| 905 | 0 | thisS = totalStuff(HI, hThick, areaT, salinity, unscale=S_unscale) |
| 906 | ||
| 907 | 0 | if (is_root_pe()) then |
| 908 | 0 | if (firstCall) then |
| 909 | 0 | totalH = thisH ; totalT = thisT ; totalS = thisS |
| 910 | 0 | write(stdout,*) 'Totals H,T,S:', thisH*HL2_unscale, thisT*T_unscale, thisS*S_unscale, ' ', mesg |
| 911 | 0 | firstCall = .false. |
| 912 | else | |
| 913 | 0 | delH = thisH - totalH |
| 914 | 0 | delT = thisT - totalT |
| 915 | 0 | delS = thisS - totalS |
| 916 | 0 | totalH = thisH ; totalT = thisT ; totalS = thisS |
| 917 | 0 | write(0,*) 'Tot/del H,T,S:', thisH*HL2_unscale, thisT*T_unscale, thisS*S_unscale, & |
| 918 | 0 | delH*HL2_unscale, delT*T_unscale, delS*S_unscale, ' ', mesg |
| 919 | endif | |
| 920 | endif | |
| 921 | ||
| 922 | 0 | end subroutine totalTandS |
| 923 | ||
| 924 | !> Returns false if the column integral of a given quantity is within roundoff | |
| 925 | 0 | logical 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 | 0 | u_sum = field(1) |
| 937 | 0 | error = 0. |
| 938 | ||
| 939 | ! Reintegrate and sum roundoff errors | |
| 940 | 0 | do k=2,nk |
| 941 | 0 | u_sum = u_sum + field(k) |
| 942 | 0 | 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 | 0 | if (present(known_answer)) then |
| 947 | 0 | expected = known_answer |
| 948 | else | |
| 949 | 0 | expected = 0. |
| 950 | endif | |
| 951 | ||
| 952 | ! Compare the column integrals against calculated roundoff error | |
| 953 | 0 | if (abs(u_sum-expected) > error) then |
| 954 | 0 | check_column_integral = .true. |
| 955 | else | |
| 956 | 0 | check_column_integral = .false. |
| 957 | endif | |
| 958 | ||
| 959 | 0 | end function check_column_integral |
| 960 | ||
| 961 | !> Returns false if the column integrals of two given quantities are within roundoff of each other | |
| 962 | 0 | logical 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 | 0 | if (present(missing_value)) then |
| 978 | 0 | misval = missing_value |
| 979 | else | |
| 980 | 0 | misval = 0. |
| 981 | endif | |
| 982 | ||
| 983 | 0 | u1_sum = field_1(1) |
| 984 | 0 | error1 = 0. |
| 985 | ||
| 986 | ! Reintegrate and sum roundoff errors | |
| 987 | 0 | do k=2,nk_1 |
| 988 | 0 | if (field_1(k) /= misval) then |
| 989 | 0 | u1_sum = u1_sum + field_1(k) |
| 990 | 0 | error1 = error1 + EPSILON(u1_sum)*MAX(ABS(u1_sum),ABS(field_1(k))) |
| 991 | endif | |
| 992 | enddo | |
| 993 | ||
| 994 | 0 | u2_sum = field_2(1) |
| 995 | 0 | error2 = 0. |
| 996 | ||
| 997 | ! Reintegrate and sum roundoff errors | |
| 998 | 0 | do k=2,nk_2 |
| 999 | 0 | if (field_2(k) /= misval) then |
| 1000 | 0 | u2_sum = u2_sum + field_2(k) |
| 1001 | 0 | 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 | 0 | if (abs(u1_sum-u2_sum) > (error1+error2)) then |
| 1007 | 0 | check_column_integrals = .true. |
| 1008 | else | |
| 1009 | 0 | check_column_integrals = .false. |
| 1010 | endif | |
| 1011 | ||
| 1012 | 0 | end function check_column_integrals |
| 1013 | ||
| 1014 | end module MOM_debugging |