diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index 465cdf2c28..515697c09e 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -2239,10 +2239,19 @@ subroutine forcing_diagnostics(fluxes, sfc_state, dt, G, diag, handles) if (handles%id_net_massout > 0 .or. handles%id_total_net_massout > 0) then do j=js,je ; do i=is,ie res(i,j) = 0.0 - if (fluxes%lprec(i,j) < 0.0) res(i,j) = res(i,j) + fluxes%lprec(i,j) - if (fluxes%vprec(i,j) < 0.0) res(i,j) = res(i,j) + fluxes%vprec(i,j) - if (fluxes%evap(i,j) < 0.0) res(i,j) = res(i,j) + fluxes%evap(i,j) - if (fluxes%seaice_melt(i,j) < 0.0) res(i,j) = res(i,j) + fluxes%seaice_melt(i,j) + if (associated(fluxes%lprec)) then + if (fluxes%lprec(i,j) < 0.0) res(i,j) = res(i,j) + fluxes%lprec(i,j) + endif + if (associated(fluxes%vprec)) then + if (fluxes%vprec(i,j) < 0.0) res(i,j) = res(i,j) + fluxes%vprec(i,j) + endif + if (associated(fluxes%evap)) then + if (fluxes%evap(i,j) < 0.0) res(i,j) = res(i,j) + fluxes%evap(i,j) + endif + if (associated(fluxes%seaice_melt)) then + if (fluxes%seaice_melt(i,j) < 0.0) & + res(i,j) = res(i,j) + fluxes%seaice_melt(i,j) + endif enddo ; enddo if (handles%id_net_massout > 0) call post_data(handles%id_net_massout, res, diag) if (handles%id_total_net_massout > 0) then @@ -2251,16 +2260,34 @@ subroutine forcing_diagnostics(fluxes, sfc_state, dt, G, diag, handles) endif endif - if (handles%id_massout_flux > 0) call post_data(handles%id_massout_flux,fluxes%netMassOut,diag) + if (handles%id_massout_flux > 0 .and. associated(fluxes%netMassOut)) & + call post_data(handles%id_massout_flux,fluxes%netMassOut,diag) if (handles%id_net_massin > 0 .or. handles%id_total_net_massin > 0) then do j=js,je ; do i=is,ie - res(i,j) = fluxes%fprec(i,j) + fluxes%lrunoff(i,j) + fluxes%frunoff(i,j) - if (fluxes%lprec(i,j) > 0.0) res(i,j) = res(i,j) + fluxes%lprec(i,j) - if (fluxes%vprec(i,j) > 0.0) res(i,j) = res(i,j) + fluxes%vprec(i,j) + res(i,j) = 0.0 + + if (associated(fluxes%fprec)) & + res(i,j) = res(i,j) + fluxes%fprec(i,j) + if (associated(fluxes%lrunoff)) & + res(i,j) = res(i,j) + fluxes%lrunoff(i,j) + if (associated(fluxes%frunoff)) & + res(i,j) = res(i,j) + fluxes%frunoff(i,j) + + if (associated(fluxes%lprec)) then + if (fluxes%lprec(i,j) > 0.0) res(i,j) = res(i,j) + fluxes%lprec(i,j) + endif + if (associated(fluxes%vprec)) then + if (fluxes%vprec(i,j) > 0.0) res(i,j) = res(i,j) + fluxes%vprec(i,j) + endif ! fluxes%cond is not needed because it is derived from %evap > 0 - if (fluxes%evap(i,j) > 0.0) res(i,j) = res(i,j) + fluxes%evap(i,j) - if (fluxes%seaice_melt(i,j) > 0.0) res(i,j) = res(i,j) + fluxes%seaice_melt(i,j) + if (associated(fluxes%evap)) then + if (fluxes%evap(i,j) > 0.0) res(i,j) = res(i,j) + fluxes%evap(i,j) + endif + if (associated(fluxes%seaice_melt)) then + if (fluxes%seaice_melt(i,j) > 0.0) & + res(i,j) = res(i,j) + fluxes%seaice_melt(i,j) + endif enddo ; enddo if (handles%id_net_massin > 0) call post_data(handles%id_net_massin, res, diag) if (handles%id_total_net_massin > 0) then @@ -2269,7 +2296,8 @@ subroutine forcing_diagnostics(fluxes, sfc_state, dt, G, diag, handles) endif endif - if (handles%id_massin_flux > 0) call post_data(handles%id_massin_flux,fluxes%netMassIn,diag) + if (handles%id_massin_flux > 0 .and. associated(fluxes%netMassIn)) & + call post_data(handles%id_massin_flux,fluxes%netMassIn,diag) if ((handles%id_evap > 0) .and. associated(fluxes%evap)) & call post_data(handles%id_evap, fluxes%evap, diag) diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index 5da7a91e17..45cfb0ac68 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -1439,7 +1439,7 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag # include "version_variable.h" character(len=40) :: mdl = "MOM_diagnostics" ! This module's name. character(len=48) :: thickness_units, flux_units - logical :: use_temperature + logical :: use_temperature, adiabatic integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, nz, nkml, nkbl integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, i, j @@ -1457,6 +1457,8 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag CS%diag => diag use_temperature = associated(tv%T) + call get_param(param_file, mdl, "ADIABATIC", adiabatic, default=.false., & + do_not_log=.true.) ! Read all relevant parameters and write them to the model log. call log_version(param_file, mdl, version) @@ -1642,10 +1644,11 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag 'Kinetic Energy Source from Horizontal Viscosity', 'm3 s-3') if (CS%id_KE_horvisc>0) call safe_alloc_ptr(CS%KE_horvisc,isd,ied,jsd,jed,nz) - CS%id_KE_dia = register_diag_field('ocean_model', 'KE_dia', diag%axesTL, Time, & - 'Kinetic Energy Source from Diapycnal Diffusion', 'm3 s-3') - if (CS%id_KE_dia>0) call safe_alloc_ptr(CS%KE_dia,isd,ied,jsd,jed,nz) - + if (.not. adiabatic) then + CS%id_KE_dia = register_diag_field('ocean_model', 'KE_dia', diag%axesTL, Time, & + 'Kinetic Energy Source from Diapycnal Diffusion', 'm3 s-3') + if (CS%id_KE_dia>0) call safe_alloc_ptr(CS%KE_dia,isd,ied,jsd,jed,nz) + endif ! gravity wave CFLs CS%id_cg1 = register_diag_field('ocean_model', 'cg1', diag%axesT1, Time, & diff --git a/src/framework/MOM_checksums.F90 b/src/framework/MOM_checksums.F90 index df014dc7a5..c6a23667db 100644 --- a/src/framework/MOM_checksums.F90 +++ b/src/framework/MOM_checksums.F90 @@ -10,11 +10,13 @@ module MOM_checksums use MOM_file_parser, only : log_version, param_file_type use MOM_hor_index, only : hor_index_type +use iso_fortran_env, only: error_unit + implicit none ; private +public :: chksum0, zchksum public :: hchksum, Bchksum, uchksum, vchksum, qchksum, is_NaN, chksum public :: hchksum_pair, uvchksum, Bchksum_pair -public :: chksum_general public :: MOM_checksums_init !> Checksums a pair of arrays (2d or 3d) staggered at tracer points @@ -72,11 +74,7 @@ module MOM_checksums module procedure is_NaN_0d, is_NaN_1d, is_NaN_2d, is_NaN_3d end interface -!> Return the bitcount of an array -interface chksum_general - module procedure chksum_general_1d, chksum_general_2d, chksum_general_3d -end interface - +integer, parameter :: bc_modulus = 1000000000 !< Modulus of checksum bitcount integer, parameter :: default_shift=0 !< The default array shift logical :: calculateStatistics=.true. !< If true, report min, max and mean. logical :: writeChksums=.true. !< If true, report the bitcount checksum @@ -85,8 +83,120 @@ module MOM_checksums contains +!> Checksum a scalar field (consistent with array checksums) +subroutine chksum0(scalar, mesg, scale, logunit) + real, intent(in) :: scalar !< The array to be checksummed + character(len=*), intent(in) :: mesg !< An identifying message + real, optional, intent(in) :: scale !< A scaling factor for this array. + integer, optional, intent(in) :: logunit !< IO unit for checksum logging + + real :: scaling !< Explicit rescaling factor + integer :: iounit !< Log IO unit + real :: rs !< Rescaled scalar + integer :: bc !< Scalar bitcount + + if (checkForNaNs .and. is_NaN(scalar)) & + call chksum_error(FATAL, 'NaN detected: '//trim(mesg)) + + scaling = 1.0 ; if (present(scale)) scaling = scale + iounit = error_unit; if(present(logunit)) iounit = logunit + + if (calculateStatistics) then + rs = scaling * scalar + if (is_root_pe()) & + call chk_sum_msg(" scalar:", rs, rs, rs, mesg, iounit) + endif + + if (.not. writeChksums) return + + bc = mod(bitcount(abs(scaling * scalar)), bc_modulus) + if (is_root_pe()) & + call chk_sum_msg(" scalar:", bc, mesg, iounit) + +end subroutine chksum0 + + +!> Checksum a 1d array (typically a column). +subroutine zchksum(array, mesg, scale, logunit) + real, dimension(:), intent(in) :: array !< The array to be checksummed + character(len=*), intent(in) :: mesg !< An identifying message + real, optional, intent(in) :: scale !< A scaling factor for this array. + integer, optional, intent(in) :: logunit !< IO unit for checksum logging + + real, allocatable, dimension(:) :: rescaled_array + real :: scaling + integer :: iounit !< Log IO unit + integer :: k + real :: aMean, aMin, aMax + integer :: bc0 + + if (checkForNaNs) then + if (is_NaN(array(:))) & + call chksum_error(FATAL, 'NaN detected: '//trim(mesg)) + endif + + scaling = 1.0 ; if (present(scale)) scaling = scale + iounit = error_unit; if(present(logunit)) iounit = logunit + + if (calculateStatistics) then + if (present(scale)) then + allocate(rescaled_array(LBOUND(array,1):UBOUND(array,1))) + rescaled_array(:) = 0.0 + do k=1, size(array, 1) + rescaled_array(k) = scale * array(k) + enddo + + call subStats(rescaled_array, aMean, aMin, aMax) + deallocate(rescaled_array) + else + call subStats(array, aMean, aMin, aMax) + endif + + if (is_root_pe()) & + call chk_sum_msg(" column:", aMean, aMin, aMax, mesg, iounit) + endif + + if (.not. writeChksums) return + + bc0 = subchk(array, scaling) + if (is_root_pe()) call chk_sum_msg(" column:", bc0, mesg, iounit) + + contains + + integer function subchk(array, scale) + real, dimension(:), intent(in) :: array !< The array to be checksummed + real, intent(in) :: scale !< A scaling factor for this array. + integer :: k, bc + subchk = 0 + do k=LBOUND(array, 1), UBOUND(array, 1) + bc = bitcount(abs(scale * array(k))) + subchk = subchk + bc + enddo + subchk=mod(subchk, bc_modulus) + end function subchk + + subroutine subStats(array, aMean, aMin, aMax) + real, dimension(:), intent(in) :: array !< The array to be checksummed + real, intent(out) :: aMean, aMin, aMax + + integer :: k, n + + aMin = array(1) + aMax = array(1) + n = 0 + do k=LBOUND(array,1), UBOUND(array,1) + aMin = min(aMin, array(k)) + aMax = max(aMax, array(k)) + n = n + 1 + enddo + aMean = sum(array(:)) / real(n) + end subroutine subStats + +end subroutine zchksum + !> Checksums on a pair of 2d arrays staggered at tracer points. -subroutine chksum_pair_h_2d(mesg, arrayA, arrayB, HI, haloshift, omit_corners, scale) +subroutine chksum_pair_h_2d(mesg, arrayA, arrayB, HI, haloshift, omit_corners, & + scale, logunit) character(len=*), intent(in) :: mesg !< Identifying messages type(hor_index_type), intent(in) :: HI !< A horizontal index type real, dimension(HI%isd:,HI%jsd:), intent(in) :: arrayA !< The first array to be checksummed @@ -94,19 +204,23 @@ subroutine chksum_pair_h_2d(mesg, arrayA, arrayB, HI, haloshift, omit_corners, s integer, optional, intent(in) :: haloshift !< The width of halos to check (default 0) logical, optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts real, optional, intent(in) :: scale !< A scaling factor for this array. + integer, optional, intent(in) :: logunit !< IO unit for checksum logging if (present(haloshift)) then - call chksum_h_2d(arrayA, 'x '//mesg, HI, haloshift, omit_corners, scale=scale) - call chksum_h_2d(arrayB, 'y '//mesg, HI, haloshift, omit_corners, scale=scale) + call chksum_h_2d(arrayA, 'x '//mesg, HI, haloshift, omit_corners, & + scale=scale, logunit=logunit) + call chksum_h_2d(arrayB, 'y '//mesg, HI, haloshift, omit_corners, & + scale=scale, logunit=logunit) else - call chksum_h_2d(arrayA, 'x '//mesg, HI, scale=scale) - call chksum_h_2d(arrayB, 'y '//mesg, HI, scale=scale) + call chksum_h_2d(arrayA, 'x '//mesg, HI, scale=scale, logunit=logunit) + call chksum_h_2d(arrayB, 'y '//mesg, HI, scale=scale, logunit=logunit) endif end subroutine chksum_pair_h_2d !> Checksums on a pair of 3d arrays staggered at tracer points. -subroutine chksum_pair_h_3d(mesg, arrayA, arrayB, HI, haloshift, omit_corners, scale) +subroutine chksum_pair_h_3d(mesg, arrayA, arrayB, HI, haloshift, omit_corners, & + scale, logunit) character(len=*), intent(in) :: mesg !< Identifying messages type(hor_index_type), intent(in) :: HI !< A horizontal index type real, dimension(HI%isd:,HI%jsd:, :), intent(in) :: arrayA !< The first array to be checksummed @@ -114,29 +228,35 @@ subroutine chksum_pair_h_3d(mesg, arrayA, arrayB, HI, haloshift, omit_corners, s integer, optional, intent(in) :: haloshift !< The width of halos to check (default 0) logical, optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts real, optional, intent(in) :: scale !< A scaling factor for this array. + integer, optional, intent(in) :: logunit !< IO unit for checksum logging if (present(haloshift)) then - call chksum_h_3d(arrayA, 'x '//mesg, HI, haloshift, omit_corners, scale=scale) - call chksum_h_3d(arrayB, 'y '//mesg, HI, haloshift, omit_corners, scale=scale) + call chksum_h_3d(arrayA, 'x '//mesg, HI, haloshift, omit_corners, & + scale=scale, logunit=logunit) + call chksum_h_3d(arrayB, 'y '//mesg, HI, haloshift, omit_corners, & + scale=scale, logunit=logunit) else - call chksum_h_3d(arrayA, 'x '//mesg, HI, scale=scale) - call chksum_h_3d(arrayB, 'y '//mesg, HI, scale=scale) + call chksum_h_3d(arrayA, 'x '//mesg, HI, scale=scale, logunit=logunit) + call chksum_h_3d(arrayB, 'y '//mesg, HI, scale=scale, logunit=logunit) endif end subroutine chksum_pair_h_3d !> Checksums a 2d array staggered at tracer points. -subroutine chksum_h_2d(array, mesg, HI, haloshift, omit_corners, scale) +subroutine chksum_h_2d(array, mesg, HI, haloshift, omit_corners, scale, logunit) type(hor_index_type), intent(in) :: HI !< A horizontal index type real, dimension(HI%isd:,HI%jsd:), intent(in) :: array !< The array to be checksummed character(len=*), intent(in) :: mesg !< An identifying message integer, optional, intent(in) :: haloshift !< The width of halos to check (default 0) logical, optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts real, optional, intent(in) :: scale !< A scaling factor for this array. + integer, optional, intent(in) :: logunit !< IO unit for checksum logging real, allocatable, dimension(:,:) :: rescaled_array real :: scaling + integer :: iounit !< Log IO unit integer :: i, j + real :: aMean, aMin, aMax integer :: bc0, bcSW, bcSE, bcNW, bcNE, hshift integer :: bcN, bcS, bcE, bcW logical :: do_corners @@ -147,20 +267,27 @@ subroutine chksum_h_2d(array, mesg, HI, haloshift, omit_corners, scale) ! if (is_NaN(array)) & ! call chksum_error(FATAL, 'NaN detected in halo: '//trim(mesg)) endif + scaling = 1.0 ; if (present(scale)) scaling = scale + iounit = error_unit; if(present(logunit)) iounit = logunit + + if (calculateStatistics) then + if (present(scale)) then + allocate( rescaled_array(LBOUND(array,1):UBOUND(array,1), & + LBOUND(array,2):UBOUND(array,2)) ) + rescaled_array(:,:) = 0.0 + do j=HI%jsc,HI%jec ; do i=HI%isc,HI%iec + rescaled_array(i,j) = scale*array(i,j) + enddo ; enddo + call subStats(HI, rescaled_array, aMean, aMin, aMax) + deallocate(rescaled_array) + else + call subStats(HI, array, aMean, aMin, aMax) + endif - if (calculateStatistics) then ; if (present(scale)) then - allocate( rescaled_array(LBOUND(array,1):UBOUND(array,1), & - LBOUND(array,2):UBOUND(array,2)) ) - rescaled_array(:,:) = 0.0 - do j=HI%jsc,HI%jec ; do i=HI%isc,HI%iec - rescaled_array(i,j) = scale*array(i,j) - enddo ; enddo - call subStats(HI, rescaled_array, mesg) - deallocate(rescaled_array) - else - call subStats(HI, array, mesg) - endif ; endif + if (is_root_pe()) & + call chk_sum_msg("h-point:", aMean, aMin, aMax, mesg, iounit) + endif if (.not.writeChksums) return @@ -179,7 +306,7 @@ subroutine chksum_h_2d(array, mesg, HI, haloshift, omit_corners, scale) bc0 = subchk(array, HI, 0, 0, scaling) if (hshift==0) then - if (is_root_pe()) call chk_sum_msg("h-point:",bc0,mesg) + if (is_root_pe()) call chk_sum_msg("h-point:", bc0, mesg, iounit) return endif @@ -191,14 +318,16 @@ subroutine chksum_h_2d(array, mesg, HI, haloshift, omit_corners, scale) bcNW = subchk(array, HI, -hshift, hshift, scaling) bcNE = subchk(array, HI, hshift, hshift, scaling) - if (is_root_pe()) call chk_sum_msg("h-point:",bc0,bcSW,bcSE,bcNW,bcNE,mesg) + if (is_root_pe()) & + call chk_sum_msg("h-point:", bc0, bcSW, bcSE, bcNW, bcNE, mesg, iounit) else bcS = subchk(array, HI, 0, -hshift, scaling) bcE = subchk(array, HI, hshift, 0, scaling) bcW = subchk(array, HI, -hshift, 0, scaling) bcN = subchk(array, HI, 0, hshift, scaling) - if (is_root_pe()) call chk_sum_msg_NSEW("h-point:",bc0,bcN,bcS,bcE,bcW,mesg) + if (is_root_pe()) & + call chk_sum_msg_NSEW("h-point:", bc0, bcN, bcS, bcE, bcW, mesg, iounit) endif contains @@ -215,16 +344,15 @@ integer function subchk(array, HI, di, dj, scale) subchk = subchk + bc enddo ; enddo call sum_across_PEs(subchk) - subchk=mod(subchk,1000000000) + subchk=mod(subchk, bc_modulus) end function subchk - subroutine subStats(HI, array, mesg) + subroutine subStats(HI, array, aMean, aMin, aMax) type(hor_index_type), intent(in) :: HI !< A horizontal index type real, dimension(HI%isd:,HI%jsd:), intent(in) :: array !< The array to be checksummed - character(len=*), intent(in) :: mesg !< An identifying message + real, intent(out) :: aMean, aMin, aMax integer :: i, j, n - real :: aMean, aMin, aMax aMin = array(HI%isc,HI%jsc) aMax = array(HI%isc,HI%jsc) @@ -239,13 +367,13 @@ subroutine subStats(HI, array, mesg) call min_across_PEs(aMin) call max_across_PEs(aMax) aMean = aMean / real(n) - if (is_root_pe()) call chk_sum_msg("h-point:",aMean,aMin,aMax,mesg) end subroutine subStats end subroutine chksum_h_2d !> Checksums on a pair of 2d arrays staggered at q-points. -subroutine chksum_pair_B_2d(mesg, arrayA, arrayB, HI, haloshift, symmetric, omit_corners, scale) +subroutine chksum_pair_B_2d(mesg, arrayA, arrayB, HI, haloshift, symmetric, & + omit_corners, scale, logunit) character(len=*), intent(in) :: mesg !< Identifying messages type(hor_index_type), intent(in) :: HI !< A horizontal index type real, dimension(HI%isd:,HI%jsd:), intent(in) :: arrayA !< The first array to be checksummed @@ -255,6 +383,7 @@ subroutine chksum_pair_B_2d(mesg, arrayA, arrayB, HI, haloshift, symmetric, omit integer, optional, intent(in) :: haloshift !< The width of halos to check (default 0) logical, optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts real, optional, intent(in) :: scale !< A scaling factor for this array. + integer, optional, intent(in) :: logunit !< IO unit for checksum logging logical :: sym @@ -262,18 +391,21 @@ subroutine chksum_pair_B_2d(mesg, arrayA, arrayB, HI, haloshift, symmetric, omit if (present(haloshift)) then call chksum_B_2d(arrayA, 'x '//mesg, HI, haloshift, symmetric=sym, & - omit_corners=omit_corners, scale=scale) + omit_corners=omit_corners, scale=scale, logunit=logunit) call chksum_B_2d(arrayB, 'y '//mesg, HI, haloshift, symmetric=sym, & - omit_corners=omit_corners, scale=scale) + omit_corners=omit_corners, scale=scale, logunit=logunit) else - call chksum_B_2d(arrayA, 'x '//mesg, HI, symmetric=sym, scale=scale) - call chksum_B_2d(arrayB, 'y '//mesg, HI, symmetric=sym, scale=scale) + call chksum_B_2d(arrayA, 'x '//mesg, HI, symmetric=sym, scale=scale, & + logunit=logunit) + call chksum_B_2d(arrayB, 'y '//mesg, HI, symmetric=sym, scale=scale, & + logunit=logunit) endif end subroutine chksum_pair_B_2d !> Checksums on a pair of 3d arrays staggered at q-points. -subroutine chksum_pair_B_3d(mesg, arrayA, arrayB, HI, haloshift, symmetric, omit_corners, scale) +subroutine chksum_pair_B_3d(mesg, arrayA, arrayB, HI, haloshift, symmetric, & + omit_corners, scale, logunit) character(len=*), intent(in) :: mesg !< Identifying messages type(hor_index_type), intent(in) :: HI !< A horizontal index type real, dimension(HI%IsdB:,HI%JsdB:, :), intent(in) :: arrayA !< The first array to be checksummed @@ -283,23 +415,27 @@ subroutine chksum_pair_B_3d(mesg, arrayA, arrayB, HI, haloshift, symmetric, omit !! symmetric computational domain. logical, optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts real, optional, intent(in) :: scale !< A scaling factor for this array. + integer, optional, intent(in) :: logunit !< IO unit for checksum logging logical :: sym if (present(haloshift)) then call chksum_B_3d(arrayA, 'x '//mesg, HI, haloshift, symmetric, & - omit_corners, scale=scale) + omit_corners, scale=scale, logunit=logunit) call chksum_B_3d(arrayB, 'y '//mesg, HI, haloshift, symmetric, & - omit_corners, scale=scale) + omit_corners, scale=scale, logunit=logunit) else - call chksum_B_3d(arrayA, 'x '//mesg, HI, symmetric=symmetric, scale=scale) - call chksum_B_3d(arrayB, 'y '//mesg, HI, symmetric=symmetric, scale=scale) + call chksum_B_3d(arrayA, 'x '//mesg, HI, symmetric=symmetric, scale=scale, & + logunit=logunit) + call chksum_B_3d(arrayB, 'y '//mesg, HI, symmetric=symmetric, scale=scale, & + logunit=logunit) endif end subroutine chksum_pair_B_3d !> Checksums a 2d array staggered at corner points. -subroutine chksum_B_2d(array, mesg, HI, haloshift, symmetric, omit_corners, scale) +subroutine chksum_B_2d(array, mesg, HI, haloshift, symmetric, omit_corners, & + scale, logunit) type(hor_index_type), intent(in) :: HI !< A horizontal index type real, dimension(HI%IsdB:,HI%JsdB:), & intent(in) :: array !< The array to be checksummed @@ -309,10 +445,13 @@ subroutine chksum_B_2d(array, mesg, HI, haloshift, symmetric, omit_corners, scal !! full symmetric computational domain. logical, optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts real, optional, intent(in) :: scale !< A scaling factor for this array. + integer, optional, intent(in) :: logunit !< IO unit for checksum logging real, allocatable, dimension(:,:) :: rescaled_array real :: scaling + integer :: iounit !< Log IO unit integer :: i, j, Is, Js + real :: aMean, aMin, aMax integer :: bc0, bcSW, bcSE, bcNW, bcNE, hshift integer :: bcN, bcS, bcE, bcW logical :: do_corners, sym, sym_stats @@ -323,24 +462,30 @@ subroutine chksum_B_2d(array, mesg, HI, haloshift, symmetric, omit_corners, scal ! if (is_NaN(array)) & ! call chksum_error(FATAL, 'NaN detected in halo: '//trim(mesg)) endif + scaling = 1.0 ; if (present(scale)) scaling = scale + iounit = error_unit; if(present(logunit)) iounit = logunit sym_stats = .false. ; if (present(symmetric)) sym_stats = symmetric if (present(haloshift)) then ; if (haloshift > 0) sym_stats = .true. ; endif - if (calculateStatistics) then ; if (present(scale)) then - allocate( rescaled_array(LBOUND(array,1):UBOUND(array,1), & - LBOUND(array,2):UBOUND(array,2)) ) - rescaled_array(:,:) = 0.0 - Is = HI%isc ; if (sym_stats) Is = HI%isc-1 - Js = HI%jsc ; if (sym_stats) Js = HI%jsc-1 - do J=Js,HI%JecB ; do I=Is,HI%IecB - rescaled_array(I,J) = scale*array(I,J) - enddo ; enddo - call subStats(HI, rescaled_array, mesg, sym_stats) - deallocate(rescaled_array) - else - call subStats(HI, array, mesg, sym_stats) - endif ; endif + if (calculateStatistics) then + if (present(scale)) then + allocate( rescaled_array(LBOUND(array,1):UBOUND(array,1), & + LBOUND(array,2):UBOUND(array,2)) ) + rescaled_array(:,:) = 0.0 + Is = HI%isc ; if (sym_stats) Is = HI%isc-1 + Js = HI%jsc ; if (sym_stats) Js = HI%jsc-1 + do J=Js,HI%JecB ; do I=Is,HI%IecB + rescaled_array(I,J) = scale*array(I,J) + enddo ; enddo + call subStats(HI, rescaled_array, sym_stats, aMean, aMin, aMax) + deallocate(rescaled_array) + else + call subStats(HI, array, sym_stats, aMean, aMin, aMax) + endif + if (is_root_pe()) & + call chk_sum_msg("B-point:", aMean, aMin, aMax, mesg, iounit) + endif if (.not.writeChksums) return @@ -361,7 +506,7 @@ subroutine chksum_B_2d(array, mesg, HI, haloshift, symmetric, omit_corners, scal sym = .false. ; if (present(symmetric)) sym = symmetric if ((hshift==0) .and. .not.sym) then - if (is_root_pe()) call chk_sum_msg("B-point:",bc0,mesg) + if (is_root_pe()) call chk_sum_msg("B-point:", bc0, mesg, iounit) return endif @@ -379,14 +524,16 @@ subroutine chksum_B_2d(array, mesg, HI, haloshift, symmetric, omit_corners, scal endif bcNE = subchk(array, HI, hshift, hshift, scaling) - if (is_root_pe()) call chk_sum_msg("B-point:",bc0,bcSW,bcSE,bcNW,bcNE,mesg) + if (is_root_pe()) & + call chk_sum_msg("B-point:", bc0, bcSW, bcSE, bcNW, bcNE, mesg, iounit) else bcS = subchk(array, HI, 0, -hshift, scaling) bcE = subchk(array, HI, hshift, 0, scaling) bcW = subchk(array, HI, -hshift, 0, scaling) bcN = subchk(array, HI, 0, hshift, scaling) - if (is_root_pe()) call chk_sum_msg_NSEW("B-point:",bc0,bcN,bcS,bcE,bcW,mesg) + if (is_root_pe()) & + call chk_sum_msg_NSEW("B-point:", bc0, bcN, bcS, bcE, bcW, mesg, iounit) endif contains @@ -405,18 +552,17 @@ integer function subchk(array, HI, di, dj, scale) subchk = subchk + bc enddo ; enddo call sum_across_PEs(subchk) - subchk=mod(subchk,1000000000) + subchk=mod(subchk, bc_modulus) end function subchk - subroutine subStats(HI, array, mesg, sym_stats) + subroutine subStats(HI, array, sym_stats, aMean, aMin, aMax) type(hor_index_type), intent(in) :: HI !< A horizontal index type real, dimension(HI%IsdB:,HI%JsdB:), intent(in) :: array !< The array to be checksummed - character(len=*), intent(in) :: mesg !< An identifying message logical, intent(in) :: sym_stats !< If true, evaluate the statistics on the !! full symmetric computational domain. + real, intent(out) :: aMean, aMin, aMax integer :: i, j, n, IsB, JsB - real :: aMean, aMin, aMax IsB = HI%isc ; if (sym_stats) IsB = HI%isc-1 JsB = HI%jsc ; if (sym_stats) JsB = HI%jsc-1 @@ -433,13 +579,13 @@ subroutine subStats(HI, array, mesg, sym_stats) call min_across_PEs(aMin) call max_across_PEs(aMax) aMean = aMean / real(n) - if (is_root_pe()) call chk_sum_msg("B-point:",aMean,aMin,aMax,mesg) end subroutine subStats end subroutine chksum_B_2d !> Checksums a pair of 2d velocity arrays staggered at C-grid locations -subroutine chksum_uv_2d(mesg, arrayU, arrayV, HI, haloshift, symmetric, omit_corners, scale) +subroutine chksum_uv_2d(mesg, arrayU, arrayV, HI, haloshift, symmetric, & + omit_corners, scale, logunit) character(len=*), intent(in) :: mesg !< Identifying messages type(hor_index_type), intent(in) :: HI !< A horizontal index type real, dimension(HI%IsdB:,HI%jsd:), intent(in) :: arrayU !< The u-component array to be checksummed @@ -449,19 +595,25 @@ subroutine chksum_uv_2d(mesg, arrayU, arrayV, HI, haloshift, symmetric, omit_cor !! symmetric computational domain. logical, optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts real, optional, intent(in) :: scale !< A scaling factor for these arrays. + integer, optional, intent(in) :: logunit !< IO unit for checksum logging if (present(haloshift)) then - call chksum_u_2d(arrayU, 'u '//mesg, HI, haloshift, symmetric, omit_corners, scale) - call chksum_v_2d(arrayV, 'v '//mesg, HI, haloshift, symmetric, omit_corners, scale) + call chksum_u_2d(arrayU, 'u '//mesg, HI, haloshift, symmetric, & + omit_corners, scale, logunit=logunit) + call chksum_v_2d(arrayV, 'v '//mesg, HI, haloshift, symmetric, & + omit_corners, scale, logunit=logunit) else - call chksum_u_2d(arrayU, 'u '//mesg, HI, symmetric=symmetric) - call chksum_v_2d(arrayV, 'v '//mesg, HI, symmetric=symmetric) + call chksum_u_2d(arrayU, 'u '//mesg, HI, symmetric=symmetric, & + logunit=logunit) + call chksum_v_2d(arrayV, 'v '//mesg, HI, symmetric=symmetric, & + logunit=logunit) endif end subroutine chksum_uv_2d !> Checksums a pair of 3d velocity arrays staggered at C-grid locations -subroutine chksum_uv_3d(mesg, arrayU, arrayV, HI, haloshift, symmetric, omit_corners, scale) +subroutine chksum_uv_3d(mesg, arrayU, arrayV, HI, haloshift, symmetric, & + omit_corners, scale, logunit) character(len=*), intent(in) :: mesg !< Identifying messages type(hor_index_type), intent(in) :: HI !< A horizontal index type real, dimension(HI%IsdB:,HI%jsd:,:), intent(in) :: arrayU !< The u-component array to be checksummed @@ -471,19 +623,25 @@ subroutine chksum_uv_3d(mesg, arrayU, arrayV, HI, haloshift, symmetric, omit_cor !! symmetric computational domain. logical, optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts real, optional, intent(in) :: scale !< A scaling factor for these arrays. + integer, optional, intent(in) :: logunit !< IO unit for checksum logging if (present(haloshift)) then - call chksum_u_3d(arrayU, 'u '//mesg, HI, haloshift, symmetric, omit_corners, scale) - call chksum_v_3d(arrayV, 'v '//mesg, HI, haloshift, symmetric, omit_corners, scale) + call chksum_u_3d(arrayU, 'u '//mesg, HI, haloshift, symmetric, & + omit_corners, scale, logunit=logunit) + call chksum_v_3d(arrayV, 'v '//mesg, HI, haloshift, symmetric, & + omit_corners, scale, logunit=logunit) else - call chksum_u_3d(arrayU, 'u '//mesg, HI, symmetric=symmetric) - call chksum_v_3d(arrayV, 'v '//mesg, HI, symmetric=symmetric) + call chksum_u_3d(arrayU, 'u '//mesg, HI, symmetric=symmetric, & + logunit=logunit) + call chksum_v_3d(arrayV, 'v '//mesg, HI, symmetric=symmetric, & + logunit=logunit) endif end subroutine chksum_uv_3d !> Checksums a 2d array staggered at C-grid u points. -subroutine chksum_u_2d(array, mesg, HI, haloshift, symmetric, omit_corners, scale) +subroutine chksum_u_2d(array, mesg, HI, haloshift, symmetric, omit_corners, & + scale, logunit) type(hor_index_type), intent(in) :: HI !< A horizontal index type real, dimension(HI%IsdB:,HI%jsd:), intent(in) :: array !< The array to be checksummed character(len=*), intent(in) :: mesg !< An identifying message @@ -492,10 +650,13 @@ subroutine chksum_u_2d(array, mesg, HI, haloshift, symmetric, omit_corners, scal !! symmetric computational domain. logical, optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts real, optional, intent(in) :: scale !< A scaling factor for this array. + integer, optional, intent(in) :: logunit !< IO unit for checksum logging real, allocatable, dimension(:,:) :: rescaled_array real :: scaling + integer :: iounit !< Log IO unit integer :: i, j, Is + real :: aMean, aMin, aMax integer :: bc0, bcSW, bcSE, bcNW, bcNE, hshift integer :: bcN, bcS, bcE, bcW logical :: do_corners, sym, sym_stats @@ -506,24 +667,30 @@ subroutine chksum_u_2d(array, mesg, HI, haloshift, symmetric, omit_corners, scal ! if (is_NaN(array)) & ! call chksum_error(FATAL, 'NaN detected in halo: '//trim(mesg)) endif - scaling = 1.0 ; if (present(scale)) scaling = scale + scaling = 1.0 ; if (present(scale)) scaling = scale + iounit = error_unit; if(present(logunit)) iounit = logunit sym_stats = .false. ; if (present(symmetric)) sym_stats = symmetric if (present(haloshift)) then ; if (haloshift > 0) sym_stats = .true. ; endif - if (calculateStatistics) then ; if (present(scale)) then - allocate( rescaled_array(LBOUND(array,1):UBOUND(array,1), & - LBOUND(array,2):UBOUND(array,2)) ) - rescaled_array(:,:) = 0.0 - Is = HI%isc ; if (sym_stats) Is = HI%isc-1 - do j=HI%jsc,HI%jec ; do I=Is,HI%IecB - rescaled_array(I,j) = scale*array(I,j) - enddo ; enddo - call subStats(HI, rescaled_array, mesg, sym_stats) - deallocate(rescaled_array) - else - call subStats(HI, array, mesg, sym_stats) - endif ; endif + if (calculateStatistics) then + if (present(scale)) then + allocate( rescaled_array(LBOUND(array,1):UBOUND(array,1), & + LBOUND(array,2):UBOUND(array,2)) ) + rescaled_array(:,:) = 0.0 + Is = HI%isc ; if (sym_stats) Is = HI%isc-1 + do j=HI%jsc,HI%jec ; do I=Is,HI%IecB + rescaled_array(I,j) = scale*array(I,j) + enddo ; enddo + call subStats(HI, rescaled_array, sym_stats, aMean, aMin, aMax) + deallocate(rescaled_array) + else + call subStats(HI, array, sym_stats, aMean, aMin, aMax) + endif + + if (is_root_pe()) & + call chk_sum_msg("u-point:", aMean, aMin, aMax, mesg, iounit) + endif if (.not.writeChksums) return @@ -544,7 +711,7 @@ subroutine chksum_u_2d(array, mesg, HI, haloshift, symmetric, omit_corners, scal sym = .false. ; if (present(symmetric)) sym = symmetric if ((hshift==0) .and. .not.sym) then - if (is_root_pe()) call chk_sum_msg("u-point:",bc0,mesg) + if (is_root_pe()) call chk_sum_msg("u-point:", bc0, mesg, iounit) return endif @@ -552,7 +719,7 @@ subroutine chksum_u_2d(array, mesg, HI, haloshift, symmetric, omit_corners, scal if (hshift==0) then bcW = subchk(array, HI, -hshift-1, 0, scaling) - if (is_root_pe()) call chk_sum_msg_W("u-point:",bc0,bcW,mesg) + if (is_root_pe()) call chk_sum_msg_W("u-point:", bc0, bcW, mesg, iounit) elseif (do_corners) then if (sym) then bcSW = subchk(array, HI, -hshift-1, -hshift, scaling) @@ -564,7 +731,8 @@ subroutine chksum_u_2d(array, mesg, HI, haloshift, symmetric, omit_corners, scal bcSE = subchk(array, HI, hshift, -hshift, scaling) bcNE = subchk(array, HI, hshift, hshift, scaling) - if (is_root_pe()) call chk_sum_msg("u-point:",bc0,bcSW,bcSE,bcNW,bcNE,mesg) + if (is_root_pe()) & + call chk_sum_msg("u-point:", bc0, bcSW, bcSE, bcNW, bcNE, mesg, iounit) else bcS = subchk(array, HI, 0, -hshift, scaling) bcE = subchk(array, HI, hshift, 0, scaling) @@ -575,7 +743,8 @@ subroutine chksum_u_2d(array, mesg, HI, haloshift, symmetric, omit_corners, scal endif bcN = subchk(array, HI, 0, hshift, scaling) - if (is_root_pe()) call chk_sum_msg_NSEW("u-point:",bc0,bcN,bcS,bcE,bcW,mesg) + if (is_root_pe()) & + call chk_sum_msg_NSEW("u-point:", bc0, bcN, bcS, bcE, bcW, mesg, iounit) endif contains @@ -594,18 +763,17 @@ integer function subchk(array, HI, di, dj, scale) subchk = subchk + bc enddo ; enddo call sum_across_PEs(subchk) - subchk=mod(subchk,1000000000) + subchk=mod(subchk, bc_modulus) end function subchk - subroutine subStats(HI, array, mesg, sym_stats) + subroutine subStats(HI, array, sym_stats, aMean, aMin, aMax) type(hor_index_type), intent(in) :: HI !< A horizontal index type real, dimension(HI%IsdB:,HI%jsd:), intent(in) :: array !< The array to be checksummed - character(len=*), intent(in) :: mesg !< An identifying message logical, intent(in) :: sym_stats !< If true, evaluate the statistics on the !! full symmetric computational domain. + real, intent(out) :: aMean, aMin, aMax integer :: i, j, n, IsB - real :: aMean, aMin, aMax IsB = HI%isc ; if (sym_stats) IsB = HI%isc-1 @@ -621,13 +789,13 @@ subroutine subStats(HI, array, mesg, sym_stats) call min_across_PEs(aMin) call max_across_PEs(aMax) aMean = aMean / real(n) - if (is_root_pe()) call chk_sum_msg("u-point:",aMean,aMin,aMax,mesg) end subroutine subStats end subroutine chksum_u_2d !> Checksums a 2d array staggered at C-grid v points. -subroutine chksum_v_2d(array, mesg, HI, haloshift, symmetric, omit_corners, scale) +subroutine chksum_v_2d(array, mesg, HI, haloshift, symmetric, omit_corners, & + scale, logunit) type(hor_index_type), intent(in) :: HI !< A horizontal index type real, dimension(HI%isd:,HI%JsdB:), intent(in) :: array !< The array to be checksummed character(len=*), intent(in) :: mesg !< An identifying message @@ -636,10 +804,13 @@ subroutine chksum_v_2d(array, mesg, HI, haloshift, symmetric, omit_corners, scal !! symmetric computational domain. logical, optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts real, optional, intent(in) :: scale !< A scaling factor for this array. + integer, optional, intent(in) :: logunit !< IO unit for checksum logging real, allocatable, dimension(:,:) :: rescaled_array real :: scaling + integer :: iounit !< Log IO unit integer :: i, j, Js + real :: aMean, aMin, aMax integer :: bc0, bcSW, bcSE, bcNW, bcNE, hshift integer :: bcN, bcS, bcE, bcW logical :: do_corners, sym, sym_stats @@ -650,24 +821,30 @@ subroutine chksum_v_2d(array, mesg, HI, haloshift, symmetric, omit_corners, scal ! if (is_NaN(array)) & ! call chksum_error(FATAL, 'NaN detected in halo: '//trim(mesg)) endif - scaling = 1.0 ; if (present(scale)) scaling = scale + scaling = 1.0 ; if (present(scale)) scaling = scale + iounit = error_unit; if(present(logunit)) iounit = logunit sym_stats = .false. ; if (present(symmetric)) sym_stats = symmetric if (present(haloshift)) then ; if (haloshift > 0) sym_stats = .true. ; endif - if (calculateStatistics) then ; if (present(scale)) then - allocate( rescaled_array(LBOUND(array,1):UBOUND(array,1), & - LBOUND(array,2):UBOUND(array,2)) ) - rescaled_array(:,:) = 0.0 - Js = HI%jsc ; if (sym_stats) Js = HI%jsc-1 - do J=Js,HI%JecB ; do i=HI%isc,HI%iec - rescaled_array(i,J) = scale*array(i,J) - enddo ; enddo - call subStats(HI, rescaled_array, mesg, sym_stats) - deallocate(rescaled_array) - else - call subStats(HI, array, mesg, sym_stats) - endif ; endif + if (calculateStatistics) then + if (present(scale)) then + allocate( rescaled_array(LBOUND(array,1):UBOUND(array,1), & + LBOUND(array,2):UBOUND(array,2)) ) + rescaled_array(:,:) = 0.0 + Js = HI%jsc ; if (sym_stats) Js = HI%jsc-1 + do J=Js,HI%JecB ; do i=HI%isc,HI%iec + rescaled_array(i,J) = scale*array(i,J) + enddo ; enddo + call subStats(HI, rescaled_array, sym_stats, aMean, aMin, aMax) + deallocate(rescaled_array) + else + call subStats(HI, array, sym_stats, aMean, aMin, aMax) + endif + + if (is_root_pe()) & + call chk_sum_msg("v-point:", aMean, aMin, aMax, mesg, iounit) + endif if (.not.writeChksums) return @@ -688,7 +865,7 @@ subroutine chksum_v_2d(array, mesg, HI, haloshift, symmetric, omit_corners, scal sym = .false. ; if (present(symmetric)) sym = symmetric if ((hshift==0) .and. .not.sym) then - if (is_root_pe()) call chk_sum_msg("v-point:",bc0,mesg) + if (is_root_pe()) call chk_sum_msg("v-point:", bc0, mesg, iounit) return endif @@ -696,7 +873,7 @@ subroutine chksum_v_2d(array, mesg, HI, haloshift, symmetric, omit_corners, scal if (hshift==0) then bcS = subchk(array, HI, 0, -hshift-1, scaling) - if (is_root_pe()) call chk_sum_msg_S("v-point:",bc0,bcS,mesg) + if (is_root_pe()) call chk_sum_msg_S("v-point:", bc0, bcS, mesg, iounit) elseif (do_corners) then if (sym) then bcSW = subchk(array, HI, -hshift, -hshift-1, scaling) @@ -708,7 +885,8 @@ subroutine chksum_v_2d(array, mesg, HI, haloshift, symmetric, omit_corners, scal bcNW = subchk(array, HI, -hshift, hshift, scaling) bcNE = subchk(array, HI, hshift, hshift, scaling) - if (is_root_pe()) call chk_sum_msg("v-point:",bc0,bcSW,bcSE,bcNW,bcNE,mesg) + if (is_root_pe()) & + call chk_sum_msg("v-point:", bc0, bcSW, bcSE, bcNW, bcNE, mesg, iounit) else if (sym) then bcS = subchk(array, HI, 0, -hshift-1, scaling) @@ -719,7 +897,8 @@ subroutine chksum_v_2d(array, mesg, HI, haloshift, symmetric, omit_corners, scal bcW = subchk(array, HI, -hshift, 0, scaling) bcN = subchk(array, HI, 0, hshift, scaling) - if (is_root_pe()) call chk_sum_msg_NSEW("v-point:",bc0,bcN,bcS,bcE,bcW,mesg) + if (is_root_pe()) & + call chk_sum_msg_NSEW("v-point:", bc0, bcN, bcS, bcE, bcW, mesg, iounit) endif contains @@ -738,18 +917,17 @@ integer function subchk(array, HI, di, dj, scale) subchk = subchk + bc enddo ; enddo call sum_across_PEs(subchk) - subchk=mod(subchk,1000000000) + subchk=mod(subchk, bc_modulus) end function subchk - subroutine subStats(HI, array, mesg, sym_stats) + subroutine subStats(HI, array, sym_stats, aMean, aMin, aMax) type(hor_index_type), intent(in) :: HI !< A horizontal index type real, dimension(HI%isd:,HI%JsdB:), intent(in) :: array !< The array to be checksummed - character(len=*), intent(in) :: mesg !< An identifying message logical, intent(in) :: sym_stats !< If true, evaluate the statistics on the !! full symmetric computational domain. + real, intent(out) :: aMean, aMin, aMax integer :: i, j, n, JsB - real :: aMean, aMin, aMax JsB = HI%jsc ; if (sym_stats) JsB = HI%jsc-1 @@ -765,23 +943,25 @@ subroutine subStats(HI, array, mesg, sym_stats) call min_across_PEs(aMin) call max_across_PEs(aMax) aMean = aMean / real(n) - if (is_root_pe()) call chk_sum_msg("v-point:",aMean,aMin,aMax,mesg) end subroutine subStats end subroutine chksum_v_2d !> Checksums a 3d array staggered at tracer points. -subroutine chksum_h_3d(array, mesg, HI, haloshift, omit_corners, scale) +subroutine chksum_h_3d(array, mesg, HI, haloshift, omit_corners, scale, logunit) type(hor_index_type), intent(in) :: HI !< A horizontal index type real, dimension(HI%isd:,HI%jsd:,:), intent(in) :: array !< The array to be checksummed character(len=*), intent(in) :: mesg !< An identifying message integer, optional, intent(in) :: haloshift !< The width of halos to check (default 0) logical, optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts real, optional, intent(in) :: scale !< A scaling factor for this array. + integer, optional, intent(in) :: logunit !< IO unit for checksum logging real, allocatable, dimension(:,:,:) :: rescaled_array real :: scaling + integer :: iounit !< Log IO unit integer :: i, j, k + real :: aMean, aMin, aMax integer :: bc0, bcSW, bcSE, bcNW, bcNE, hshift integer :: bcN, bcS, bcE, bcW logical :: do_corners @@ -792,22 +972,29 @@ subroutine chksum_h_3d(array, mesg, HI, haloshift, omit_corners, scale) ! if (is_NaN(array)) & ! call chksum_error(FATAL, 'NaN detected in halo: '//trim(mesg)) endif - scaling = 1.0 ; if (present(scale)) scaling = scale - if (calculateStatistics) then ; if (present(scale)) then - allocate( rescaled_array(LBOUND(array,1):UBOUND(array,1), & - LBOUND(array,2):UBOUND(array,2), & - LBOUND(array,3):UBOUND(array,3)) ) - rescaled_array(:,:,:) = 0.0 - do k=1,size(array,3) ; do j=HI%jsc,HI%jec ; do i=HI%isc,HI%iec - rescaled_array(i,j,k) = scale*array(i,j,k) - enddo ; enddo ; enddo + scaling = 1.0 ; if (present(scale)) scaling = scale + iounit = error_unit; if(present(logunit)) iounit = logunit + + if (calculateStatistics) then + if (present(scale)) then + allocate( rescaled_array(LBOUND(array,1):UBOUND(array,1), & + LBOUND(array,2):UBOUND(array,2), & + LBOUND(array,3):UBOUND(array,3)) ) + rescaled_array(:,:,:) = 0.0 + do k=1,size(array,3) ; do j=HI%jsc,HI%jec ; do i=HI%isc,HI%iec + rescaled_array(i,j,k) = scale*array(i,j,k) + enddo ; enddo ; enddo + + call subStats(HI, rescaled_array, aMean, aMin, aMax) + deallocate(rescaled_array) + else + call subStats(HI, array, aMean, aMin, aMax) + endif - call subStats(HI, rescaled_array, mesg) - deallocate(rescaled_array) - else - call subStats(HI, array, mesg) - endif ; endif + if (is_root_pe()) & + call chk_sum_msg("h-point:", aMean, aMin, aMax, mesg, iounit) + endif if (.not.writeChksums) return @@ -826,7 +1013,7 @@ subroutine chksum_h_3d(array, mesg, HI, haloshift, omit_corners, scale) bc0 = subchk(array, HI, 0, 0, scaling) if (hshift==0) then - if (is_root_pe()) call chk_sum_msg("h-point:",bc0,mesg) + if (is_root_pe()) call chk_sum_msg("h-point:", bc0, mesg, iounit) return endif @@ -838,14 +1025,16 @@ subroutine chksum_h_3d(array, mesg, HI, haloshift, omit_corners, scale) bcNW = subchk(array, HI, -hshift, hshift, scaling) bcNE = subchk(array, HI, hshift, hshift, scaling) - if (is_root_pe()) call chk_sum_msg("h-point:",bc0,bcSW,bcSE,bcNW,bcNE,mesg) + if (is_root_pe()) & + call chk_sum_msg("h-point:", bc0, bcSW, bcSE, bcNW, bcNE, mesg, iounit) else bcS = subchk(array, HI, 0, -hshift, scaling) bcE = subchk(array, HI, hshift, 0, scaling) bcW = subchk(array, HI, -hshift, 0, scaling) bcN = subchk(array, HI, 0, hshift, scaling) - if (is_root_pe()) call chk_sum_msg_NSEW("h-point:",bc0,bcN,bcS,bcE,bcW,mesg) + if (is_root_pe()) & + call chk_sum_msg_NSEW("h-point:", bc0, bcN, bcS, bcE, bcW, mesg, iounit) endif contains @@ -863,16 +1052,15 @@ integer function subchk(array, HI, di, dj, scale) subchk = subchk + bc enddo ; enddo ; enddo call sum_across_PEs(subchk) - subchk=mod(subchk,1000000000) + subchk=mod(subchk, bc_modulus) end function subchk - subroutine subStats(HI, array, mesg) + subroutine subStats(HI, array, aMean, aMin, aMax) type(hor_index_type), intent(in) :: HI !< A horizontal index type real, dimension(HI%isd:,HI%jsd:,:), intent(in) :: array !< The array to be checksummed - character(len=*), intent(in) :: mesg !< An identifying message + real, intent(out) :: aMean, aMin, aMax integer :: i, j, k, n - real :: aMean, aMin, aMax aMin = array(HI%isc,HI%jsc,1) aMax = array(HI%isc,HI%jsc,1) @@ -887,13 +1075,13 @@ subroutine subStats(HI, array, mesg) call min_across_PEs(aMin) call max_across_PEs(aMax) aMean = aMean / real(n) - if (is_root_pe()) call chk_sum_msg("h-point:",aMean,aMin,aMax,mesg) end subroutine subStats end subroutine chksum_h_3d !> Checksums a 3d array staggered at corner points. -subroutine chksum_B_3d(array, mesg, HI, haloshift, symmetric, omit_corners, scale) +subroutine chksum_B_3d(array, mesg, HI, haloshift, symmetric, omit_corners, & + scale, logunit) type(hor_index_type), intent(in) :: HI !< A horizontal index type real, dimension(HI%IsdB:,HI%JsdB:,:), intent(in) :: array !< The array to be checksummed character(len=*), intent(in) :: mesg !< An identifying message @@ -902,10 +1090,13 @@ subroutine chksum_B_3d(array, mesg, HI, haloshift, symmetric, omit_corners, scal !! symmetric computational domain. logical, optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts real, optional, intent(in) :: scale !< A scaling factor for this array. + integer, optional, intent(in) :: logunit !< IO unit for checksum logging real, allocatable, dimension(:,:,:) :: rescaled_array real :: scaling + integer :: iounit !< Log IO unit integer :: i, j, k, Is, Js + real :: aMean, aMin, aMax integer :: bc0, bcSW, bcSE, bcNW, bcNE, hshift integer :: bcN, bcS, bcE, bcW logical :: do_corners, sym, sym_stats @@ -916,25 +1107,32 @@ subroutine chksum_B_3d(array, mesg, HI, haloshift, symmetric, omit_corners, scal ! if (is_NaN(array)) & ! call chksum_error(FATAL, 'NaN detected in halo: '//trim(mesg)) endif + scaling = 1.0 ; if (present(scale)) scaling = scale + iounit = error_unit; if(present(logunit)) iounit = logunit sym_stats = .false. ; if (present(symmetric)) sym_stats = symmetric if (present(haloshift)) then ; if (haloshift > 0) sym_stats = .true. ; endif - if (calculateStatistics) then ; if (present(scale)) then - allocate( rescaled_array(LBOUND(array,1):UBOUND(array,1), & - LBOUND(array,2):UBOUND(array,2), & - LBOUND(array,3):UBOUND(array,3)) ) - rescaled_array(:,:,:) = 0.0 - Is = HI%isc ; if (sym_stats) Is = HI%isc-1 - Js = HI%jsc ; if (sym_stats) Js = HI%jsc-1 - do k=1,size(array,3) ; do J=Js,HI%JecB ; do I=Is,HI%IecB - rescaled_array(I,J,k) = scale*array(I,J,k) - enddo ; enddo ; enddo - call subStats(HI, rescaled_array, mesg, sym_stats) - deallocate(rescaled_array) - else - call subStats(HI, array, mesg, sym_stats) - endif ; endif + if (calculateStatistics) then + if (present(scale)) then + allocate( rescaled_array(LBOUND(array,1):UBOUND(array,1), & + LBOUND(array,2):UBOUND(array,2), & + LBOUND(array,3):UBOUND(array,3)) ) + rescaled_array(:,:,:) = 0.0 + Is = HI%isc ; if (sym_stats) Is = HI%isc-1 + Js = HI%jsc ; if (sym_stats) Js = HI%jsc-1 + do k=1,size(array,3) ; do J=Js,HI%JecB ; do I=Is,HI%IecB + rescaled_array(I,J,k) = scale*array(I,J,k) + enddo ; enddo ; enddo + call subStats(HI, rescaled_array, sym_stats, aMean, aMin, aMax) + deallocate(rescaled_array) + else + call subStats(HI, array, sym_stats, aMean, aMin, aMax) + endif + + if (is_root_pe()) & + call chk_sum_msg("B-point:", aMean, aMin, aMax, mesg, iounit) + endif if (.not.writeChksums) return @@ -955,7 +1153,7 @@ subroutine chksum_B_3d(array, mesg, HI, haloshift, symmetric, omit_corners, scal sym = .false. ; if (present(symmetric)) sym = symmetric if ((hshift==0) .and. .not.sym) then - if (is_root_pe()) call chk_sum_msg("B-point:",bc0,mesg) + if (is_root_pe()) call chk_sum_msg("B-point:", bc0, mesg, iounit) return endif @@ -973,7 +1171,8 @@ subroutine chksum_B_3d(array, mesg, HI, haloshift, symmetric, omit_corners, scal endif bcNE = subchk(array, HI, hshift, hshift, scaling) - if (is_root_pe()) call chk_sum_msg("B-point:",bc0,bcSW,bcSE,bcNW,bcNE,mesg) + if (is_root_pe()) & + call chk_sum_msg("B-point:", bc0, bcSW, bcSE, bcNW, bcNE, mesg, iounit) else if (sym) then bcS = subchk(array, HI, 0, -hshift-1, scaling) @@ -985,7 +1184,8 @@ subroutine chksum_B_3d(array, mesg, HI, haloshift, symmetric, omit_corners, scal bcE = subchk(array, HI, hshift, 0, scaling) bcN = subchk(array, HI, 0, hshift, scaling) - if (is_root_pe()) call chk_sum_msg_NSEW("B-point:",bc0,bcN,bcS,bcE,bcW,mesg) + if (is_root_pe()) & + call chk_sum_msg_NSEW("B-point:", bc0, bcN, bcS, bcE, bcW, mesg, iounit) endif contains @@ -1004,18 +1204,17 @@ integer function subchk(array, HI, di, dj, scale) subchk = subchk + bc enddo ; enddo ; enddo call sum_across_PEs(subchk) - subchk=mod(subchk,1000000000) + subchk=mod(subchk, bc_modulus) end function subchk - subroutine subStats(HI, array, mesg, sym_stats) + subroutine subStats(HI, array, sym_stats, aMean, aMin, aMax) type(hor_index_type), intent(in) :: HI !< A horizontal index type real, dimension(HI%IsdB:,HI%JsdB:,:), intent(in) :: array !< The array to be checksummed - character(len=*), intent(in) :: mesg !< An identifying message logical, intent(in) :: sym_stats !< If true, evaluate the statistics on the !! full symmetric computational domain. + real, intent(out) :: aMean, aMin, aMax integer :: i, j, k, n, IsB, JsB - real :: aMean, aMin, aMax IsB = HI%isc ; if (sym_stats) IsB = HI%isc-1 JsB = HI%jsc ; if (sym_stats) JsB = HI%jsc-1 @@ -1031,13 +1230,13 @@ subroutine subStats(HI, array, mesg, sym_stats) call min_across_PEs(aMin) call max_across_PEs(aMax) aMean = aMean / real(n) - if (is_root_pe()) call chk_sum_msg("B-point:",aMean,aMin,aMax,mesg) end subroutine subStats end subroutine chksum_B_3d !> Checksums a 3d array staggered at C-grid u points. -subroutine chksum_u_3d(array, mesg, HI, haloshift, symmetric, omit_corners, scale) +subroutine chksum_u_3d(array, mesg, HI, haloshift, symmetric, omit_corners, & + scale, logunit) type(hor_index_type), intent(in) :: HI !< A horizontal index type real, dimension(HI%isdB:,HI%Jsd:,:), intent(in) :: array !< The array to be checksummed character(len=*), intent(in) :: mesg !< An identifying message @@ -1046,10 +1245,13 @@ subroutine chksum_u_3d(array, mesg, HI, haloshift, symmetric, omit_corners, scal !! symmetric computational domain. logical, optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts real, optional, intent(in) :: scale !< A scaling factor for this array. + integer, optional, intent(in) :: logunit !< IO unit for checksum logging real, allocatable, dimension(:,:,:) :: rescaled_array real :: scaling + integer :: iounit !< Log IO unit integer :: i, j, k, Is + real :: aMean, aMin, aMax integer :: bc0, bcSW, bcSE, bcNW, bcNE, hshift integer :: bcN, bcS, bcE, bcW logical :: do_corners, sym, sym_stats @@ -1060,24 +1262,30 @@ subroutine chksum_u_3d(array, mesg, HI, haloshift, symmetric, omit_corners, scal ! if (is_NaN(array)) & ! call chksum_error(FATAL, 'NaN detected in halo: '//trim(mesg)) endif + scaling = 1.0 ; if (present(scale)) scaling = scale + iounit = error_unit; if(present(logunit)) iounit = logunit sym_stats = .false. ; if (present(symmetric)) sym_stats = symmetric if (present(haloshift)) then ; if (haloshift > 0) sym_stats = .true. ; endif - if (calculateStatistics) then ; if (present(scale)) then - allocate( rescaled_array(LBOUND(array,1):UBOUND(array,1), & - LBOUND(array,2):UBOUND(array,2), & - LBOUND(array,3):UBOUND(array,3)) ) - rescaled_array(:,:,:) = 0.0 - Is = HI%isc ; if (sym_stats) Is = HI%isc-1 - do k=1,size(array,3) ; do j=HI%jsc,HI%jec ; do I=Is,HI%IecB - rescaled_array(I,j,k) = scale*array(I,j,k) - enddo ; enddo ; enddo - call subStats(HI, rescaled_array, mesg, sym_stats) - deallocate(rescaled_array) - else - call subStats(HI, array, mesg, sym_stats) - endif ; endif + if (calculateStatistics) then + if (present(scale)) then + allocate( rescaled_array(LBOUND(array,1):UBOUND(array,1), & + LBOUND(array,2):UBOUND(array,2), & + LBOUND(array,3):UBOUND(array,3)) ) + rescaled_array(:,:,:) = 0.0 + Is = HI%isc ; if (sym_stats) Is = HI%isc-1 + do k=1,size(array,3) ; do j=HI%jsc,HI%jec ; do I=Is,HI%IecB + rescaled_array(I,j,k) = scale*array(I,j,k) + enddo ; enddo ; enddo + call subStats(HI, rescaled_array, sym_stats, aMean, aMin, aMax) + deallocate(rescaled_array) + else + call subStats(HI, array, sym_stats, aMean, aMin, aMax) + endif + if (is_root_pe()) & + call chk_sum_msg("u-point:", aMean, aMin, aMax, mesg, iounit) + endif if (.not.writeChksums) return @@ -1098,7 +1306,7 @@ subroutine chksum_u_3d(array, mesg, HI, haloshift, symmetric, omit_corners, scal sym = .false. ; if (present(symmetric)) sym = symmetric if ((hshift==0) .and. .not.sym) then - if (is_root_pe()) call chk_sum_msg("u-point:",bc0,mesg) + if (is_root_pe()) call chk_sum_msg("u-point:", bc0, mesg, iounit) return endif @@ -1106,7 +1314,7 @@ subroutine chksum_u_3d(array, mesg, HI, haloshift, symmetric, omit_corners, scal if (hshift==0) then bcW = subchk(array, HI, -hshift-1, 0, scaling) - if (is_root_pe()) call chk_sum_msg_W("u-point:",bc0,bcW,mesg) + if (is_root_pe()) call chk_sum_msg_W("u-point:", bc0, bcW, mesg, iounit) elseif (do_corners) then if (sym) then bcSW = subchk(array, HI, -hshift-1, -hshift, scaling) @@ -1118,7 +1326,8 @@ subroutine chksum_u_3d(array, mesg, HI, haloshift, symmetric, omit_corners, scal bcSE = subchk(array, HI, hshift, -hshift, scaling) bcNE = subchk(array, HI, hshift, hshift, scaling) - if (is_root_pe()) call chk_sum_msg("u-point:",bc0,bcSW,bcSE,bcNW,bcNE,mesg) + if (is_root_pe()) & + call chk_sum_msg("u-point:", bc0, bcSW, bcSE, bcNW, bcNE, mesg, iounit) else bcS = subchk(array, HI, 0, -hshift, scaling) bcE = subchk(array, HI, hshift, 0, scaling) @@ -1129,7 +1338,8 @@ subroutine chksum_u_3d(array, mesg, HI, haloshift, symmetric, omit_corners, scal endif bcN = subchk(array, HI, 0, hshift, scaling) - if (is_root_pe()) call chk_sum_msg_NSEW("u-point:",bc0,bcN,bcS,bcE,bcW,mesg) + if (is_root_pe()) & + call chk_sum_msg_NSEW("u-point:", bc0, bcN, bcS, bcE, bcW, mesg, iounit) endif contains @@ -1148,18 +1358,17 @@ integer function subchk(array, HI, di, dj, scale) subchk = subchk + bc enddo ; enddo ; enddo call sum_across_PEs(subchk) - subchk=mod(subchk,1000000000) + subchk=mod(subchk, bc_modulus) end function subchk - subroutine subStats(HI, array, mesg, sym_stats) + subroutine subStats(HI, array, sym_stats, aMean, aMin, aMax) type(hor_index_type), intent(in) :: HI !< A horizontal index type real, dimension(HI%IsdB:,HI%jsd:,:), intent(in) :: array !< The array to be checksummed - character(len=*), intent(in) :: mesg !< An identifying message logical, intent(in) :: sym_stats !< If true, evaluate the statistics on the !! full symmetric computational domain. + real, intent(out) :: aMean, aMin, aMax integer :: i, j, k, n, IsB - real :: aMean, aMin, aMax IsB = HI%isc ; if (sym_stats) IsB = HI%isc-1 @@ -1175,85 +1384,13 @@ subroutine subStats(HI, array, mesg, sym_stats) call min_across_PEs(aMin) call max_across_PEs(aMax) aMean = aMean / real(n) - if (is_root_pe()) call chk_sum_msg("u-point:",aMean,aMin,aMax,mesg) end subroutine subStats end subroutine chksum_u_3d -!> Return the bitcount of an arbitrarily sized 3d array -integer function chksum_general_3d( array, scale_factor, istart, iend, jstart, jend, kstart, kend ) & - result(subchk) - real, dimension(:,:,:), intent(in) :: array !< Array to be checksummed - real, optional, intent(in) :: scale_factor !< Factor to scale array by before checksum - integer, optional, intent(in) :: istart !< Starting index in the i-direction - integer, optional, intent(in) :: iend !< Ending index in the i-direction - integer, optional, intent(in) :: jstart !< Starting index in the j-direction - integer, optional, intent(in) :: jend !< Ending index in the j-direction - integer, optional, intent(in) :: kstart !< Starting index in the k-direction - integer, optional, intent(in) :: kend !< Ending index in the k-direction - integer :: i, j, k, bc, is, ie, js, je, ks, ke - real :: scale - - ! By default do not scale - scale = 1. - if (present(scale_factor)) scale = scale_factor - - ! Set the loop indices based on full array - is = LBOUND(array,1) ; ie = UBOUND(array,1) - js = LBOUND(array,2) ; je = UBOUND(array,2) - ks = LBOUND(array,3) ; ke = UBOUND(array,3) - - ! Override indices if subdomain requested - if (present(istart)) is = istart ; if (present(iend)) ie = iend - if (present(jstart)) js = jstart ; if (present(jend)) je = jend - if (present(kstart)) ks = kstart ; if (present(kend)) ke = kend - - subchk = 0 - do k=ks,ke ; do j=js,je ; do i=is,ie - bc = bitcount(abs(scale*array(i,j,k))) - subchk = subchk + bc - enddo ; enddo ; enddo - call sum_across_PEs(subchk) - subchk=mod(subchk,1000000000) -end function chksum_general_3d - -!> Return the bitcount of an arbitrarily sized 2d array by promotion to a 3d array -integer function chksum_general_2d( array_2d, scale_factor, istart, iend, jstart, jend ) - real, dimension(:,:), intent(in) :: array_2d !< Array to be checksummed - real, optional, intent(in) :: scale_factor !< Factor to scale array by before checksum - integer, optional, intent(in) :: istart !< Starting index in the i-direction - integer, optional, intent(in) :: iend !< Ending index in the i-direction - integer, optional, intent(in) :: jstart !< Starting index in the j-direction - integer, optional, intent(in) :: jend !< Ending index in the j-direction - integer :: is, ie, js, je - real, dimension(:,:,:), allocatable :: array_3d !< Promotion from 2d to 3d array - - is = LBOUND(array_2d,1) ; ie = UBOUND(array_2d,1) - js = LBOUND(array_2d,2) ; je = UBOUND(array_2d,2) - allocate(array_3d(is:ie, js:je,1)) - array_3d(:,:,1) = array_2d(:,:) - chksum_general_2d = chksum_general_3d( array_3d, scale_factor, istart, iend, jstart, jend ) - deallocate(array_3d) -end function chksum_general_2d - -!> Return the bitcount of an arbitrarily sized 1d array by promotion to a 3d array -integer function chksum_general_1d( array_1d, scale_factor, istart, iend ) - real, dimension(:), intent(in) :: array_1d !< Array to be checksummed - real, optional, intent(in) :: scale_factor !< Factor to scale array by before checksum - integer, optional, intent(in) :: istart !< Starting index in the i-direction - integer, optional, intent(in) :: iend !< Ending index in the i-direction - integer :: is, ie - real, dimension(:,:,:), allocatable :: array_3d !< Promotion from 2d to 3d array - - is = LBOUND(array_1d,1) ; ie = UBOUND(array_1d,1) - allocate(array_3d(is:ie, 1,1)) - array_3d(:,1,1) = array_1d(:) - chksum_general_1d = chksum_general_3d( array_3d, scale_factor, istart, iend ) - deallocate(array_3d) -end function chksum_general_1d - !> Checksums a 3d array staggered at C-grid v points. -subroutine chksum_v_3d(array, mesg, HI, haloshift, symmetric, omit_corners, scale) +subroutine chksum_v_3d(array, mesg, HI, haloshift, symmetric, omit_corners, & + scale, logunit) type(hor_index_type), intent(in) :: HI !< A horizontal index type real, dimension(HI%isd:,HI%JsdB:,:), intent(in) :: array !< The array to be checksummed character(len=*), intent(in) :: mesg !< An identifying message @@ -1262,12 +1399,15 @@ subroutine chksum_v_3d(array, mesg, HI, haloshift, symmetric, omit_corners, scal !! symmetric computational domain. logical, optional, intent(in) :: omit_corners !< If true, avoid checking diagonal shifts real, optional, intent(in) :: scale !< A scaling factor for this array. + integer, optional, intent(in) :: logunit !< IO unit for checksum logging real, allocatable, dimension(:,:,:) :: rescaled_array real :: scaling + integer :: iounit !< Log IO unit integer :: i, j, k, Js integer :: bc0, bcSW, bcSE, bcNW, bcNE, hshift integer :: bcN, bcS, bcE, bcW + real :: aMean, aMin, aMax logical :: do_corners, sym, sym_stats if (checkForNaNs) then @@ -1276,24 +1416,30 @@ subroutine chksum_v_3d(array, mesg, HI, haloshift, symmetric, omit_corners, scal ! if (is_NaN(array)) & ! call chksum_error(FATAL, 'NaN detected in halo: '//trim(mesg)) endif + scaling = 1.0 ; if (present(scale)) scaling = scale + iounit = error_unit; if(present(logunit)) iounit = logunit sym_stats = .false. ; if (present(symmetric)) sym_stats = symmetric if (present(haloshift)) then ; if (haloshift > 0) sym_stats = .true. ; endif - if (calculateStatistics) then ; if (present(scale)) then - allocate( rescaled_array(LBOUND(array,1):UBOUND(array,1), & - LBOUND(array,2):UBOUND(array,2), & - LBOUND(array,3):UBOUND(array,3)) ) - rescaled_array(:,:,:) = 0.0 - Js = HI%jsc ; if (sym_stats) Js = HI%jsc-1 - do k=1,size(array,3) ; do J=Js,HI%JecB ; do i=HI%isc,HI%iec - rescaled_array(i,J,k) = scale*array(i,J,k) - enddo ; enddo ; enddo - call subStats(HI, rescaled_array, mesg, sym_stats) - deallocate(rescaled_array) - else - call subStats(HI, array, mesg, sym_stats) - endif ; endif + if (calculateStatistics) then + if (present(scale)) then + allocate( rescaled_array(LBOUND(array,1):UBOUND(array,1), & + LBOUND(array,2):UBOUND(array,2), & + LBOUND(array,3):UBOUND(array,3)) ) + rescaled_array(:,:,:) = 0.0 + Js = HI%jsc ; if (sym_stats) Js = HI%jsc-1 + do k=1,size(array,3) ; do J=Js,HI%JecB ; do i=HI%isc,HI%iec + rescaled_array(i,J,k) = scale*array(i,J,k) + enddo ; enddo ; enddo + call subStats(HI, rescaled_array, sym_stats, aMean, aMin, aMax) + deallocate(rescaled_array) + else + call subStats(HI, array, sym_stats, aMean, aMin, aMax) + endif + if (is_root_pe()) & + call chk_sum_msg("v-point:", aMean, aMin, aMax, mesg, iounit) + endif if (.not.writeChksums) return @@ -1314,7 +1460,7 @@ subroutine chksum_v_3d(array, mesg, HI, haloshift, symmetric, omit_corners, scal sym = .false. ; if (present(symmetric)) sym = symmetric if ((hshift==0) .and. .not.sym) then - if (is_root_pe()) call chk_sum_msg("v-point:",bc0,mesg) + if (is_root_pe()) call chk_sum_msg("v-point:", bc0, mesg, iounit) return endif @@ -1322,7 +1468,7 @@ subroutine chksum_v_3d(array, mesg, HI, haloshift, symmetric, omit_corners, scal if (hshift==0) then bcS = subchk(array, HI, 0, -hshift-1, scaling) - if (is_root_pe()) call chk_sum_msg_S("v-point:",bc0,bcS,mesg) + if (is_root_pe()) call chk_sum_msg_S("v-point:", bc0, bcS, mesg, iounit) elseif (do_corners) then if (sym) then bcSW = subchk(array, HI, -hshift, -hshift-1, scaling) @@ -1334,7 +1480,8 @@ subroutine chksum_v_3d(array, mesg, HI, haloshift, symmetric, omit_corners, scal bcNW = subchk(array, HI, -hshift, hshift, scaling) bcNE = subchk(array, HI, hshift, hshift, scaling) - if (is_root_pe()) call chk_sum_msg("v-point:",bc0,bcSW,bcSE,bcNW,bcNE,mesg) + if (is_root_pe()) & + call chk_sum_msg("v-point:", bc0, bcSW, bcSE, bcNW, bcNE, mesg, iounit) else if (sym) then bcS = subchk(array, HI, 0, -hshift-1, scaling) @@ -1345,7 +1492,8 @@ subroutine chksum_v_3d(array, mesg, HI, haloshift, symmetric, omit_corners, scal bcW = subchk(array, HI, -hshift, 0, scaling) bcN = subchk(array, HI, 0, hshift, scaling) - if (is_root_pe()) call chk_sum_msg_NSEW("v-point:",bc0,bcN,bcS,bcE,bcW,mesg) + if (is_root_pe()) & + call chk_sum_msg_NSEW("v-point:", bc0, bcN, bcS, bcE, bcW, mesg, iounit) endif contains @@ -1364,18 +1512,18 @@ integer function subchk(array, HI, di, dj, scale) subchk = subchk + bc enddo ; enddo ; enddo call sum_across_PEs(subchk) - subchk=mod(subchk,1000000000) + subchk=mod(subchk, bc_modulus) end function subchk - subroutine subStats(HI, array, mesg, sym_stats) + !subroutine subStats(HI, array, mesg, sym_stats) + subroutine subStats(HI, array, sym_stats, aMean, aMin, aMax) type(hor_index_type), intent(in) :: HI !< A horizontal index type real, dimension(HI%isd:,HI%JsdB:,:), intent(in) :: array !< The array to be checksummed - character(len=*), intent(in) :: mesg !< An identifying message logical, intent(in) :: sym_stats !< If true, evaluate the statistics on the !! full symmetric computational domain. + real, intent(out) :: aMean, aMin, aMax !< Mean/min/max of array over domain integer :: i, j, k, n, JsB - real :: aMean, aMin, aMax JsB = HI%jsc ; if (sym_stats) JsB = HI%jsc-1 @@ -1391,7 +1539,6 @@ subroutine subStats(HI, array, mesg, sym_stats) call min_across_PEs(aMin) call max_across_PEs(aMax) aMean = aMean / real(n) - if (is_root_pe()) call chk_sum_msg("v-point:",aMean,aMin,aMax,mesg) end subroutine subStats end subroutine chksum_v_3d @@ -1590,15 +1737,18 @@ function is_NaN_3d(x) end function is_NaN_3d !> Write a message including the checksum of the non-shifted array -subroutine chk_sum_msg1(fmsg,bc0,mesg) +subroutine chk_sum_msg1(fmsg, bc0, mesg, iounit) character(len=*), intent(in) :: fmsg !< A checksum code-location specific preamble character(len=*), intent(in) :: mesg !< An identifying message supplied by top-level caller integer, intent(in) :: bc0 !< The bitcount of the non-shifted array - if (is_root_pe()) write(0,'(A,1(A,I10,X),A)') fmsg," c=",bc0,trim(mesg) + integer, intent(in) :: iounit !< Checksum logger IO unit + + if (is_root_pe()) & + write(iounit, '(A,1(A,I10,X),A)') fmsg, " c=", bc0, trim(mesg) end subroutine chk_sum_msg1 !> Write a message including checksums of non-shifted and diagonally shifted arrays -subroutine chk_sum_msg5(fmsg,bc0,bcSW,bcSE,bcNW,bcNE,mesg) +subroutine chk_sum_msg5(fmsg, bc0, bcSW, bcSE, bcNW, bcNE, mesg, iounit) character(len=*), intent(in) :: fmsg !< A checksum code-location specific preamble character(len=*), intent(in) :: mesg !< An identifying message supplied by top-level caller integer, intent(in) :: bc0 !< The bitcount of the non-shifted array @@ -1606,12 +1756,14 @@ subroutine chk_sum_msg5(fmsg,bc0,bcSW,bcSE,bcNW,bcNE,mesg) integer, intent(in) :: bcSE !< The bitcount for SE shifted array integer, intent(in) :: bcNW !< The bitcount for NW shifted array integer, intent(in) :: bcNE !< The bitcount for NE shifted array - if (is_root_pe()) write(0,'(A,5(A,I10,1X),A)') & - fmsg," c=",bc0,"sw=",bcSW,"se=",bcSE,"nw=",bcNW,"ne=",bcNE,trim(mesg) + integer, intent(in) :: iounit !< Checksum logger IO unit + + if (is_root_pe()) write(iounit, '(A,5(A,I10,1X),A)') & + fmsg, " c=", bc0, "sw=", bcSW, "se=", bcSE, "nw=", bcNW, "ne=", bcNE, trim(mesg) end subroutine chk_sum_msg5 !> Write a message including checksums of non-shifted and laterally shifted arrays -subroutine chk_sum_msg_NSEW(fmsg,bc0,bcN,bcS,bcE,bcW,mesg) +subroutine chk_sum_msg_NSEW(fmsg, bc0, bcN, bcS, bcE, bcW, mesg, iounit) character(len=*), intent(in) :: fmsg !< A checksum code-location specific preamble character(len=*), intent(in) :: mesg !< An identifying message supplied by top-level caller integer, intent(in) :: bc0 !< The bitcount of the non-shifted array @@ -1619,49 +1771,59 @@ subroutine chk_sum_msg_NSEW(fmsg,bc0,bcN,bcS,bcE,bcW,mesg) integer, intent(in) :: bcS !< The bitcount for S shifted array integer, intent(in) :: bcE !< The bitcount for E shifted array integer, intent(in) :: bcW !< The bitcount for W shifted array - if (is_root_pe()) write(0,'(A,5(A,I10,1X),A)') & - fmsg," c=",bc0,"N=",bcN,"S=",bcS,"E=",bcE,"W=",bcW,trim(mesg) + integer, intent(in) :: iounit !< Checksum logger IO unit + + if (is_root_pe()) write(iounit, '(A,5(A,I10,1X),A)') & + fmsg, " c=", bc0, "N=", bcN, "S=", bcS, "E=", bcE, "W=", bcW, trim(mesg) end subroutine chk_sum_msg_NSEW !> Write a message including checksums of non-shifted and southward shifted arrays -subroutine chk_sum_msg_S(fmsg,bc0,bcS,mesg) +subroutine chk_sum_msg_S(fmsg, bc0, bcS, mesg, iounit) character(len=*), intent(in) :: fmsg !< A checksum code-location specific preamble character(len=*), intent(in) :: mesg !< An identifying message supplied by top-level caller integer, intent(in) :: bc0 !< The bitcount of the non-shifted array integer, intent(in) :: bcS !< The bitcount of the south-shifted array - if (is_root_pe()) write(0,'(A,2(A,I10,1X),A)') & - fmsg," c=",bc0,"S=",bcS,trim(mesg) + integer, intent(in) :: iounit !< Checksum logger IO unit + + if (is_root_pe()) write(iounit, '(A,2(A,I10,1X),A)') & + fmsg, " c=", bc0, "S=", bcS, trim(mesg) end subroutine chk_sum_msg_S !> Write a message including checksums of non-shifted and westward shifted arrays -subroutine chk_sum_msg_W(fmsg,bc0,bcW,mesg) +subroutine chk_sum_msg_W(fmsg, bc0, bcW, mesg, iounit) character(len=*), intent(in) :: fmsg !< A checksum code-location specific preamble character(len=*), intent(in) :: mesg !< An identifying message supplied by top-level caller integer, intent(in) :: bc0 !< The bitcount of the non-shifted array integer, intent(in) :: bcW !< The bitcount of the west-shifted array - if (is_root_pe()) write(0,'(A,2(A,I10,1X),A)') & - fmsg," c=",bc0,"W=",bcW,trim(mesg) + integer, intent(in) :: iounit !< Checksum logger IO unit + + if (is_root_pe()) write(iounit, '(A,2(A,I10,1X),A)') & + fmsg, " c=", bc0, "W=", bcW, trim(mesg) end subroutine chk_sum_msg_W !> Write a message including checksums of non-shifted and southwestward shifted arrays -subroutine chk_sum_msg2(fmsg,bc0,bcSW,mesg) +subroutine chk_sum_msg2(fmsg, bc0, bcSW, mesg, iounit) character(len=*), intent(in) :: fmsg !< A checksum code-location specific preamble character(len=*), intent(in) :: mesg !< An identifying message supplied by top-level caller integer, intent(in) :: bc0 !< The bitcount of the non-shifted array integer, intent(in) :: bcSW !< The bitcount of the southwest-shifted array - if (is_root_pe()) write(0,'(A,2(A,I9,1X),A)') & - fmsg," c=",bc0,"s/w=",bcSW,trim(mesg) + integer, intent(in) :: iounit !< Checksum logger IO unit + + if (is_root_pe()) write(iounit, '(A,2(A,I9,1X),A)') & + fmsg, " c=", bc0, "s/w=", bcSW, trim(mesg) end subroutine chk_sum_msg2 !> Write a message including the global mean, maximum and minimum of an array -subroutine chk_sum_msg3(fmsg,aMean,aMin,aMax,mesg) +subroutine chk_sum_msg3(fmsg, aMean, aMin, aMax, mesg, iounit) character(len=*), intent(in) :: fmsg !< A checksum code-location specific preamble character(len=*), intent(in) :: mesg !< An identifying message supplied by top-level caller real, intent(in) :: aMean !< The mean value of the array real, intent(in) :: aMin !< The minimum value of the array real, intent(in) :: aMax !< The maximum value of the array - if (is_root_pe()) write(0,'(A,3(A,ES25.16,1X),A)') & - fmsg," mean=",aMean,"min=",aMin,"max=",aMax,trim(mesg) + integer, intent(in) :: iounit !< Checksum logger IO unit + + if (is_root_pe()) write(iounit, '(A,3(A,ES25.16,1X),A)') & + fmsg, " mean=", aMean, "min=", aMin, "max=", aMax, trim(mesg) end subroutine chk_sum_msg3 !> MOM_checksums_init initializes the MOM_checksums module. As it happens, the @@ -1686,7 +1848,7 @@ end subroutine chksum_error !> Does a bitcount of a number by first casting to an integer and then using BTEST !! to check bit by bit -integer function bitcount( x ) +integer function bitcount(x) real :: x !< Number to be bitcount ! Local variables diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index 5acc240445..9320f503b5 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -4,7 +4,8 @@ module MOM_diag_mediator ! This file is part of MOM6. See LICENSE.md for the license. -use MOM_checksums, only : chksum_general +use MOM_checksums, only : chksum0, zchksum +use MOM_checksums, only : hchksum, uchksum, vchksum, Bchksum use MOM_coms, only : PE_here use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end use MOM_cpu_clock, only : CLOCK_MODULE, CLOCK_ROUTINE @@ -12,6 +13,7 @@ module MOM_diag_mediator use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_grid, only : ocean_grid_type use MOM_io, only : slasher, vardesc, query_vardesc, mom_read_data +use MOM_io, only : get_filename_appendix use MOM_safe_alloc, only : safe_alloc_ptr, safe_alloc_alloc use MOM_string_functions, only : lowercase use MOM_time_manager, only : time_type @@ -237,7 +239,7 @@ module MOM_diag_mediator type, public :: diag_ctrl integer :: available_diag_doc_unit = -1 !< The unit number of a diagnostic documentation file. !! This file is open if available_diag_doc_unit is > 0. - integer :: chksum_diag_doc_unit = -1 !< The unit number of a diagnostic documentation file. + integer :: chksum_iounit = -1 !< The unit number of a diagnostic documentation file. !! This file is open if available_diag_doc_unit is > 0. logical :: diag_as_chksum !< If true, log chksums in a text file instead of posting diagnostics @@ -327,6 +329,9 @@ module MOM_diag_mediator real, dimension(:,:,:), allocatable :: h_old #endif + !> Number of checksum-only diagnostics + integer :: num_chksum_diags + end type diag_ctrl ! CPU clocks @@ -1210,7 +1215,9 @@ subroutine post_data_0d(diag_field_id, field, diag_cs, is_static) 'post_data_0d: Unregistered diagnostic id') diag => diag_cs%diags(diag_field_id) do while (associated(diag)) - if (is_stat) then + if (diag_cs%diag_as_chksum) then + call chksum0(field, diag%debug_str, logunit=diag_cs%chksum_iounit) + else if (is_stat) then used = send_data(diag%fms_diag_id, field) elseif (diag_cs%ave_enabled) then used = send_data(diag%fms_diag_id, field, diag_cs%time_end) @@ -1260,7 +1267,9 @@ subroutine post_data_1d_k(diag_field_id, field, diag_cs, is_static) locfield => field endif - if (is_stat) then + if (diag_cs%diag_as_chksum) then + call zchksum(locfield, diag%debug_str, logunit=diag_cs%chksum_iounit) + else if (is_stat) then used = send_data(diag%fms_diag_id, locfield) elseif (diag_cs%ave_enabled) then used = send_data(diag%fms_diag_id, locfield, diag_cs%time_end, weight=diag_cs%time_int) @@ -1397,9 +1406,20 @@ subroutine post_data_2d_low(diag, field, diag_cs, is_static, mask) endif if (diag_cs%diag_as_chksum) then - chksum = chksum_general(locfield) - if (is_root_pe()) then - call log_chksum_diag(diag_cs%chksum_diag_doc_unit, diag%debug_str, chksum) + if (diag%axes%is_h_point) then + call hchksum(locfield, diag%debug_str, diag_cs%G%HI, & + logunit=diag_cs%chksum_iounit) + else if (diag%axes%is_u_point) then + call uchksum(locfield, diag%debug_str, diag_cs%G%HI, & + logunit=diag_cs%chksum_iounit) + else if (diag%axes%is_v_point) then + call vchksum(locfield, diag%debug_str, diag_cs%G%HI, & + logunit=diag_cs%chksum_iounit) + else if (diag%axes%is_q_point) then + call Bchksum(locfield, diag%debug_str, diag_cs%G%HI, & + logunit=diag_cs%chksum_iounit) + else + call MOM_error(FATAL, "post_data_2d_low: unknown axis type.") endif else if (is_stat) then @@ -1672,9 +1692,20 @@ subroutine post_data_3d_low(diag, field, diag_cs, is_static, mask) if (diag%fms_diag_id>0) then if (diag_cs%diag_as_chksum) then - chksum = chksum_general(locfield) - if (is_root_pe()) then - call log_chksum_diag(diag_cs%chksum_diag_doc_unit, diag%debug_str, chksum) + if (diag%axes%is_h_point) then + call hchksum(locfield, diag%debug_str, diag_cs%G%HI, & + logunit=diag_cs%chksum_iounit) + else if (diag%axes%is_u_point) then + call uchksum(locfield, diag%debug_str, diag_cs%G%HI, & + logunit=diag_cs%chksum_iounit) + else if (diag%axes%is_v_point) then + call vchksum(locfield, diag%debug_str, diag_cs%G%HI, & + logunit=diag_cs%chksum_iounit) + else if (diag%axes%is_q_point) then + call Bchksum(locfield, diag%debug_str, diag_cs%G%HI, & + logunit=diag_cs%chksum_iounit) + else + call MOM_error(FATAL, "post_data_3d_low: unknown axis type.") endif else if (is_stat) then @@ -1759,8 +1790,13 @@ subroutine post_xy_average(diag_cs, diag, field) averaged_field, averaged_mask) endif - used = send_data(diag%fms_xyave_diag_id, averaged_field, diag_cs%time_end, & - weight=diag_cs%time_int, mask=averaged_mask) + if (diag_cs%diag_as_chksum) then + call zchksum(averaged_field, trim(diag%debug_str)//'_xyave', & + logunit=diag_CS%chksum_iounit) + else + used = send_data(diag%fms_xyave_diag_id, averaged_field, diag_cs%time_end, & + weight=diag_cs%time_int, mask=averaged_mask) + endif end subroutine post_xy_average !> This subroutine enables the accumulation of time averages over the specified time interval. @@ -1951,6 +1987,9 @@ integer function register_diag_field(module_name, field_name, axes_in, init_time !Register downsampled diagnostics do dl=2,MAX_DSAMP_LEV + ! Do not attempt to checksum the downsampled diagnostics + if (diag_cs%diag_as_chksum) cycle + new_module_name = trim(module_name)//'_d2' if (axes_in%rank == 3 .or. axes_in%rank == 2 ) then @@ -2115,9 +2154,10 @@ logical function register_diag_field_expand_cmor(dm_id, module_name, field_name, range=range, mask_variant=mask_variant, standard_name=standard_name, & verbose=verbose, do_not_log=do_not_log, err_msg=err_msg, & interp_method=interp_method, tile_count=tile_count) - call attach_cell_methods(fms_id, axes, cm_string, & - cell_methods, x_cell_method, y_cell_method, v_cell_method, & - v_extensive=v_extensive) + if (.not. diag_cs%diag_as_chksum) & + call attach_cell_methods(fms_id, axes, cm_string, cell_methods, & + x_cell_method, y_cell_method, v_cell_method, & + v_extensive=v_extensive) if (is_root_pe() .and. diag_CS%available_diag_doc_unit > 0) then msg = '' if (present(cmor_field_name)) msg = 'CMOR equivalent is "'//trim(cmor_field_name)//'"' @@ -2133,8 +2173,9 @@ logical function register_diag_field_expand_cmor(dm_id, module_name, field_name, range=range, mask_variant=mask_variant, standard_name=standard_name, & verbose=verbose, do_not_log=do_not_log, err_msg=err_msg, & interp_method=interp_method, tile_count=tile_count) - call attach_cell_methods(fms_xyave_id, axes%xyave_axes, cm_string, & - cell_methods, v_cell_method, v_extensive=v_extensive) + if (.not. diag_cs%diag_as_chksum) & + call attach_cell_methods(fms_xyave_id, axes%xyave_axes, cm_string, & + cell_methods, v_cell_method, v_extensive=v_extensive) if (is_root_pe() .and. diag_CS%available_diag_doc_unit > 0) then msg = '' if (present(cmor_field_name)) msg = 'CMOR equivalent is "'//trim(cmor_field_name)//'_xyave"' @@ -2154,7 +2195,7 @@ logical function register_diag_field_expand_cmor(dm_id, module_name, field_name, endif ! For the CMOR variation of the above diagnostic - if (present(cmor_field_name)) then + if (present(cmor_field_name) .and. .not. diag_cs%diag_as_chksum) then ! Fallback values for strings set to "NULL" posted_cmor_units = "not provided" ! posted_cmor_standard_name = "not provided" ! Values might be able to be replaced with a CS%missing field? @@ -2251,7 +2292,10 @@ integer function register_diag_field_expand_axes(module_name, field_name, axes, volume_id = axes%id_volume ! Get the FMS diagnostic id - if (present(interp_method) .or. axes%is_h_point) then + if (axes%diag_cs%diag_as_chksum) then + fms_id = axes%diag_cs%num_chksum_diags + 1 + axes%diag_cs%num_chksum_diags = fms_id + else if (present(interp_method) .or. axes%is_h_point) then ! If interp_method is provided we must use it if (area_id>0) then if (volume_id>0) then @@ -2563,9 +2607,16 @@ function register_scalar_field(module_name, field_name, init_time, diag_cs, & diag => null() cmor_diag => null() - fms_id = register_diag_field_fms(module_name, field_name, init_time, & - long_name=long_name, units=units, missing_value=MOM_missing_value, & - range=range, standard_name=standard_name, do_not_log=do_not_log, err_msg=err_msg) + if (diag_cs%diag_as_chksum) then + fms_id = diag_cs%num_chksum_diags + 1 + diag_cs%num_chksum_diags = fms_id + else + fms_id = register_diag_field_fms(module_name, field_name, init_time, & + long_name=long_name, units=units, missing_value=MOM_missing_value, & + range=range, standard_name=standard_name, do_not_log=do_not_log, & + err_msg=err_msg) + endif + if (fms_id /= DIAG_FIELD_NOT_FOUND) then dm_id = get_new_diag_id(diag_cs) call alloc_diag_with_id(dm_id, diag_cs, diag) @@ -2669,11 +2720,17 @@ function register_static_field(module_name, field_name, axes, & diag => null() cmor_diag => null() - fms_id = register_static_field_fms(module_name, field_name, axes%handles, & - long_name=long_name, units=units, missing_value=MOM_missing_value, & - range=range, mask_variant=mask_variant, standard_name=standard_name, & - do_not_log=do_not_log, & - interp_method=interp_method, tile_count=tile_count, area=area) + if (diag_cs%diag_as_chksum) then + fms_id = diag_cs%num_chksum_diags + 1 + diag_cs%num_chksum_diags = fms_id + else + fms_id = register_static_field_fms(module_name, field_name, axes%handles, & + long_name=long_name, units=units, missing_value=MOM_missing_value, & + range=range, mask_variant=mask_variant, standard_name=standard_name, & + do_not_log=do_not_log, & + interp_method=interp_method, tile_count=tile_count, area=area) + endif + if (fms_id /= DIAG_FIELD_NOT_FOUND) then dm_id = get_new_diag_id(diag_cs) call alloc_diag_with_id(dm_id, diag_cs, diag) @@ -2681,20 +2738,28 @@ function register_static_field(module_name, field_name, axes, & diag%fms_diag_id = fms_id diag%debug_str = trim(module_name)//"-"//trim(field_name) if (present(conversion)) diag%conversion_factor = conversion - if (present(x_cell_method)) then - call get_diag_axis_name(axes%handles(1), axis_name) - call diag_field_add_attribute(fms_id, 'cell_methods', trim(axis_name)//':'//trim(x_cell_method)) - endif - if (present(y_cell_method)) then - call get_diag_axis_name(axes%handles(2), axis_name) - call diag_field_add_attribute(fms_id, 'cell_methods', trim(axis_name)//':'//trim(y_cell_method)) - endif - if (present(area_cell_method)) then - call diag_field_add_attribute(fms_id, 'cell_methods', 'area:'//trim(area_cell_method)) + + if (diag_cs%diag_as_chksum) then + diag%axes => axes + else + if (present(x_cell_method)) then + call get_diag_axis_name(axes%handles(1), axis_name) + call diag_field_add_attribute(fms_id, 'cell_methods', & + trim(axis_name)//':'//trim(x_cell_method)) + endif + if (present(y_cell_method)) then + call get_diag_axis_name(axes%handles(2), axis_name) + call diag_field_add_attribute(fms_id, 'cell_methods', & + trim(axis_name)//':'//trim(y_cell_method)) + endif + if (present(area_cell_method)) then + call diag_field_add_attribute(fms_id, 'cell_methods', & + 'area:'//trim(area_cell_method)) + endif endif endif - if (present(cmor_field_name)) then + if (present(cmor_field_name) .and. .not. diag_cs%diag_as_chksum) then ! Fallback values for strings set to "not provided" posted_cmor_units = "not provided" posted_cmor_standard_name = "not provided" @@ -2905,6 +2970,7 @@ subroutine diag_mediator_init(G, GV, US, nz, param_file, diag_cs, doc_file_dir) ! This include declares and sets the variable "version". #include "version_variable.h" character(len=40) :: mdl = "MOM_diag_mediator" ! This module's name. + character(len=32) :: filename_appendix = '' !fms appendix to filename for ensemble runs id_clock_diag_mediator = cpu_clock_id('(Ocean diagnostics framework)', grain=CLOCK_MODULE) id_clock_diag_remap = cpu_clock_id('(Ocean diagnostics remapping)', grain=CLOCK_ROUTINE) @@ -2955,6 +3021,9 @@ subroutine diag_mediator_init(G, GV, US, nz, param_file, diag_cs, doc_file_dir) 'a text file containing the checksum (bitcount) of the array.', & default=.false.) + if (diag_cs%diag_as_chksum) & + diag_cs%num_chksum_diags = 0 + ! Keep pointers grid, h, T, S needed diagnostic remapping diag_cs%G => G diag_cs%GV => GV @@ -3023,15 +3092,25 @@ subroutine diag_mediator_init(G, GV, US, nz, param_file, diag_cs, doc_file_dir) endif endif - if (is_root_pe() .and. (diag_CS%chksum_diag_doc_unit < 0) .and. diag_CS%diag_as_chksum) then - write(this_pe,'(i6.6)') PE_here() - doc_file_dflt = "chksum_diag."//this_pe + if (is_root_pe() .and. (diag_CS%chksum_iounit < 0) .and. diag_CS%diag_as_chksum) then + !write(this_pe,'(i6.6)') PE_here() + !doc_file_dflt = "chksum_diag."//this_pe + doc_file_dflt = "chksum_diag" call get_param(param_file, mdl, "CHKSUM_DIAG_FILE", doc_file, & "A file into which to write all checksums of the "//& "diagnostics listed in the diag_table.", & - default=doc_file_dflt, do_not_log=(diag_CS%chksum_diag_doc_unit/=-1)) + default=doc_file_dflt, do_not_log=(diag_CS%chksum_iounit/=-1)) + + call get_filename_appendix(filename_appendix) + if (len_trim(filename_appendix) > 0) then + doc_file = trim(doc_file) //'.'//trim(filename_appendix) + endif +#ifdef STATSLABEL + doc_file = trim(doc_file)//"."//trim(adjustl(STATSLABEL)) +#endif + if (len_trim(doc_file) > 0) then - new_file = .true. ; if (diag_CS%chksum_diag_doc_unit /= -1) new_file = .false. + new_file = .true. ; if (diag_CS%chksum_iounit /= -1) new_file = .false. ! Find an unused unit number. do new_unit=512,42,-1 inquire( new_unit, opened=opened) @@ -3045,16 +3124,16 @@ subroutine diag_mediator_init(G, GV, US, nz, param_file, diag_cs, doc_file_dir) doc_path = trim(slasher(doc_file_dir))//trim(doc_file) endif ; endif - diag_CS%chksum_diag_doc_unit = new_unit + diag_CS%chksum_iounit = new_unit if (new_file) then - open(diag_CS%chksum_diag_doc_unit, file=trim(doc_path), access='SEQUENTIAL', form='FORMATTED', & + open(diag_CS%chksum_iounit, file=trim(doc_path), access='SEQUENTIAL', form='FORMATTED', & action='WRITE', status='REPLACE', iostat=ios) else ! This file is being reopened, and should be appended. - open(diag_CS%chksum_diag_doc_unit, file=trim(doc_path), access='SEQUENTIAL', form='FORMATTED', & + open(diag_CS%chksum_iounit, file=trim(doc_path), access='SEQUENTIAL', form='FORMATTED', & action='WRITE', status='OLD', position='APPEND', iostat=ios) endif - inquire(diag_CS%chksum_diag_doc_unit, opened=opened) + inquire(diag_CS%chksum_iounit, opened=opened) if ((.not.opened) .or. (ios /= 0)) then call MOM_error(FATAL, "Failed to open checksum diags file "//trim(doc_path)//".") endif @@ -3205,8 +3284,8 @@ subroutine diag_mediator_end(time, diag_CS, end_diag_manager) if (diag_CS%available_diag_doc_unit > -1) then close(diag_CS%available_diag_doc_unit) ; diag_CS%available_diag_doc_unit = -3 endif - if (diag_CS%chksum_diag_doc_unit > -1) then - close(diag_CS%chksum_diag_doc_unit) ; diag_CS%chksum_diag_doc_unit = -3 + if (diag_CS%chksum_iounit > -1) then + close(diag_CS%chksum_iounit) ; diag_CS%chksum_iounit = -3 endif deallocate(diag_cs%diags) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index e21192bfae..25d4eadb7d 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -2938,32 +2938,34 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di enddo endif - CS%id_Tdif = register_diag_field('ocean_model',"Tflx_dia_diff",diag%axesTi, & - Time, "Diffusive diapycnal temperature flux across interfaces", & - "degC m s-1") - CS%id_Tadv = register_diag_field('ocean_model',"Tflx_dia_adv",diag%axesTi, & - Time, "Advective diapycnal temperature flux across interfaces", & - "degC m s-1") - CS%id_Sdif = register_diag_field('ocean_model',"Sflx_dia_diff",diag%axesTi, & - Time, "Diffusive diapycnal salnity flux across interfaces", & - "psu m s-1") - CS%id_Sadv = register_diag_field('ocean_model',"Sflx_dia_adv",diag%axesTi, & - Time, "Advective diapycnal salnity flux across interfaces", & - "psu m s-1") - CS%id_MLD_003 = register_diag_field('ocean_model', 'MLD_003', diag%axesT1, Time, & - 'Mixed layer depth (delta rho = 0.03)', 'm', conversion=US%Z_to_m, & - cmor_field_name='mlotst', cmor_long_name='Ocean Mixed Layer Thickness Defined by Sigma T', & - cmor_standard_name='ocean_mixed_layer_thickness_defined_by_sigma_t') - CS%id_mlotstsq = register_diag_field('ocean_model','mlotstsq',diag%axesT1, Time, & - long_name='Square of Ocean Mixed Layer Thickness Defined by Sigma T', & - standard_name='square_of_ocean_mixed_layer_thickness_defined_by_sigma_t', & - units='m2', conversion=US%Z_to_m**2) - CS%id_MLD_0125 = register_diag_field('ocean_model','MLD_0125',diag%axesT1,Time, & - 'Mixed layer depth (delta rho = 0.125)', 'm', conversion=US%Z_to_m) - CS%id_subMLN2 = register_diag_field('ocean_model','subML_N2',diag%axesT1,Time, & - 'Squared buoyancy frequency below mixed layer', 's-2') - CS%id_MLD_user = register_diag_field('ocean_model','MLD_user',diag%axesT1,Time, & - 'Mixed layer depth (used defined)', 'm', conversion=US%Z_to_m) + if (use_temperature) then + CS%id_Tdif = register_diag_field('ocean_model',"Tflx_dia_diff",diag%axesTi, & + Time, "Diffusive diapycnal temperature flux across interfaces", & + "degC m s-1") + CS%id_Tadv = register_diag_field('ocean_model',"Tflx_dia_adv",diag%axesTi, & + Time, "Advective diapycnal temperature flux across interfaces", & + "degC m s-1") + CS%id_Sdif = register_diag_field('ocean_model',"Sflx_dia_diff",diag%axesTi, & + Time, "Diffusive diapycnal salnity flux across interfaces", & + "psu m s-1") + CS%id_Sadv = register_diag_field('ocean_model',"Sflx_dia_adv",diag%axesTi, & + Time, "Advective diapycnal salnity flux across interfaces", & + "psu m s-1") + CS%id_MLD_003 = register_diag_field('ocean_model', 'MLD_003', diag%axesT1, Time, & + 'Mixed layer depth (delta rho = 0.03)', 'm', conversion=US%Z_to_m, & + cmor_field_name='mlotst', cmor_long_name='Ocean Mixed Layer Thickness Defined by Sigma T', & + cmor_standard_name='ocean_mixed_layer_thickness_defined_by_sigma_t') + CS%id_mlotstsq = register_diag_field('ocean_model','mlotstsq',diag%axesT1, Time, & + long_name='Square of Ocean Mixed Layer Thickness Defined by Sigma T', & + standard_name='square_of_ocean_mixed_layer_thickness_defined_by_sigma_t', & + units='m2', conversion=US%Z_to_m**2) + CS%id_MLD_0125 = register_diag_field('ocean_model','MLD_0125',diag%axesT1,Time, & + 'Mixed layer depth (delta rho = 0.125)', 'm', conversion=US%Z_to_m) + CS%id_subMLN2 = register_diag_field('ocean_model','subML_N2',diag%axesT1,Time, & + 'Squared buoyancy frequency below mixed layer', 's-2') + CS%id_MLD_user = register_diag_field('ocean_model','MLD_user',diag%axesT1,Time, & + 'Mixed layer depth (used defined)', 'm', conversion=US%Z_to_m) + endif call get_param(param_file, mdl, "DIAG_MLD_DENSITY_DIFF", CS%MLDdensityDifference, & "The density difference used to determine a diagnostic mixed "//& "layer depth, MLD_user, following the definition of Levitus 1982. "//&