MOM_checksum_packages.F90

Go to the documentation of this file.

module mom_checksum_packages

!***********************************************************************
!*                   GNU General Public License                        *
!* This file is a part of MOM.                                         *
!*                                                                     *
!* MOM is free software; you can redistribute it and/or modify it and  *
!* are expected to follow the terms of the GNU General Public License  *
!* as published by the Free Software Foundation; either version 2 of   *
!* the License, or (at your option) any later version.                 *
!*                                                                     *
!* MOM is distributed in the hope that it will be useful, but WITHOUT  *
!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY  *
!* or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public    *
!* License for more details.                                           *
!*                                                                     *
!* For the full text of the GNU General Public License,                *
!* write to: Free Software Foundation, Inc.,                           *
!*           675 Mass Ave, Cambridge, MA 02139, USA.                   *
!* or see:   http://www.gnu.org/licenses/gpl.html                      *
!***********************************************************************

!   This module provdes a several routines that do check-sums of groups
! of variables in the various dynamic solver routines.

use mom_debugging, only : hchksum, uvchksum
use mom_domains, only : sum_across_pes, min_across_pes, max_across_pes
use mom_error_handler, only : mom_mesg, is_root_pe
use mom_grid, only : ocean_grid_type
use mom_variables, only : thermo_var_ptrs
use mom_verticalgrid, only : verticalgrid_type

implicit none ; private

public mom_state_chksum, mom_thermo_chksum, mom_accel_chksum
public mom_state_stats

interface mom_state_chksum
  module procedure mom_state_chksum_5arg
  module procedure mom_state_chksum_3arg
end interface

#include <MOM_memory.h>

type :: stats
  private
  real :: minimum = 1.e34, maximum = -1.e34, average = 0.
end type stats

contains

! =============================================================================

subroutine mom_state_chksum_5arg(mesg, u, v, h, uh, vh, G, GV, haloshift, symmetric)
  character(len=*),                          &
                           intent(in) :: mesg !< A message that appears on the chksum lines.
  type(ocean_grid_type),   intent(in) :: G    !< The ocean's grid structure.
  type(verticalgrid_type), intent(in) :: GV   !< The ocean's vertical grid structure.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                           intent(in) :: u    !< The zonal velocity, in m s-1.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                           intent(in) :: v    !< The meridional velocity, in m s-1.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  &
                           intent(in) :: h    !< Layer thicknesses, in H (usually m or kg m-2).
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                           intent(in) :: uh   !< Volume flux through zonal faces = u*h*dy, m3 s-1.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                           intent(in) :: vh   !< Volume flux through meridional
                                              !! faces = v*h*dx, in m3 s-1.
  integer, optional,       intent(in) :: haloshift
  logical, optional,       intent(in) :: symmetric
!   This subroutine writes out chksums for the model's basic state variables.
! Arguments: mesg - A message that appears on the chksum lines.
!  (in)      u - Zonal velocity, in m s-1.
!  (in)      v - Meridional velocity, in m s-1.
!  (in)      h - Layer thickness, in m.
!  (in)      uh - Volume flux through zonal faces = u*h*dy, m3 s-1.
!  (in)      vh - Volume flux through meridional faces = v*h*dx, in m3 s-1.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
  integer :: is, ie, js, je, nz, hs
  logical :: sym
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke

  ! Note that for the chksum calls to be useful for reproducing across PE
  ! counts, there must be no redundant points, so all variables use is..ie
  ! and js...je as their extent.
  hs=1; if (present(haloshift)) hs=haloshift
  sym=.false.; if (present(symmetric)) sym=symmetric
  call uvchksum(mesg//" [uv]", u, v, g%HI, haloshift=hs, symmetric=sym)
  call hchksum(h, mesg//" h", g%HI, haloshift=hs, scale=gv%H_to_m)
  call uvchksum(mesg//" [uv]h", uh, vh, g%HI, haloshift=hs, &
                symmetric=sym, scale=gv%H_to_m)
end subroutine mom_state_chksum_5arg

! =============================================================================

subroutine mom_state_chksum_3arg(mesg, u, v, h, G, GV, haloshift, symmetric)
  character(len=*),        intent(in) :: mesg !< A message that appears on the chksum lines.
  type(ocean_grid_type),   intent(in) :: G    !< The ocean's grid structure.
  type(verticalgrid_type), intent(in) :: GV   !< The ocean's vertical grid structure.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                           intent(in) :: u    !< Zonal velocity, in m s-1.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                           intent(in) :: v    !< Meridional velocity, in m s-1.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  &
                           intent(in) :: h    !< Layer thicknesses, in H (usually m or kg m-2).
  integer, optional,       intent(in) :: haloshift
  logical, optional,       intent(in) :: symmetric
!   This subroutine writes out chksums for the model's basic state variables.
! Arguments: mesg - A message that appears on the chksum lines.
!  (in)      u - Zonal velocity, in m s-1.
!  (in)      v - Meridional velocity, in m s-1.
!  (in)      h - Layer thickness, in m.
!  (in)      uh - Volume flux through zonal faces = u*h*dy, m3 s-1.
!  (in)      vh - Volume flux through meridional faces = v*h*dx, in m3 s-1.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
  integer :: is, ie, js, je, nz, hs
  logical :: sym
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke

  ! Note that for the chksum calls to be useful for reproducing across PE
  ! counts, there must be no redundant points, so all variables use is..ie
  ! and js...je as their extent.
  hs=1; if (present(haloshift)) hs=haloshift
  sym=.false.; if (present(symmetric)) sym=symmetric
  call uvchksum(mesg//" u", u, v, g%HI,haloshift=hs, symmetric=sym)
  call hchksum(h, mesg//" h",g%HI, haloshift=hs, scale=gv%H_to_m)
end subroutine mom_state_chksum_3arg

! =============================================================================

subroutine mom_thermo_chksum(mesg, tv, G, haloshift)
  character(len=*),         intent(in) :: mesg !< A message that appears on the chksum lines.
  type(thermo_var_ptrs),    intent(in) :: tv   !< A structure pointing to various
                                               !! thermodynamic variables.
  type(ocean_grid_type),    intent(in) :: G    !< The ocean's grid structure.
  integer, optional,        intent(in) :: haloshift
!   This subroutine writes out chksums for the model's thermodynamic state
! variables.
! Arguments: mesg - A message that appears on the chksum lines.
!  (in)      tv - A structure containing pointers to any thermodynamic
!                 fields that are in use.
!  (in)      G - The ocean's grid structure.
  integer :: is, ie, js, je, nz, hs
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  hs=1; if (present(haloshift)) hs=haloshift

  if (associated(tv%T)) call hchksum(tv%T, mesg//" T",g%HI,haloshift=hs)
  if (associated(tv%S)) call hchksum(tv%S, mesg//" S",g%HI,haloshift=hs)
  if (associated(tv%frazil)) call hchksum(tv%frazil, mesg//" frazil",g%HI,haloshift=hs)
  if (associated(tv%salt_deficit)) call hchksum(tv%salt_deficit, mesg//" salt deficit",g%HI,haloshift=hs)

end subroutine mom_thermo_chksum

! =============================================================================

subroutine mom_accel_chksum(mesg, CAu, CAv, PFu, PFv, diffu, diffv, G, GV, pbce, &
                            u_accel_bt, v_accel_bt, symmetric)
  character(len=*),         intent(in) :: mesg !< A message that appears on the chksum lines.
  type(ocean_grid_type),    intent(in) :: G    !< The ocean's grid structure.
  type(verticalgrid_type),  intent(in) :: GV   !< The ocean's vertical grid structure.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                            intent(in) :: CAu  !< Zonal acceleration due to Coriolis
                                               !! and momentum advection terms, in m s-2.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                            intent(in) :: CAv  !< Meridional acceleration due to Coriolis
                                               !! and momentum advection terms, in m s-2.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                            intent(in) :: PFu  !< Zonal acceleration due to pressure gradients
                                               !! (equal to -dM/dx) in m s-2.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                            intent(in) :: PFv  !< Meridional acceleration due to pressure gradients
                                               !! (equal to -dM/dy) in m s-2.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                            intent(in) :: diffu !< Zonal acceleration due to convergence of the
                                                !! along-isopycnal stress tensor, in m s-2.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                            intent(in) :: diffv !< Meridional acceleration due to convergence of
                                                !! the along-isopycnal stress tensor, in m s-2.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  &
                  optional, intent(in) :: pbce !< The baroclinic pressure anomaly in each layer
                                               !! due to free surface height anomalies, in
                                               !! m2 s-2 H-1.                                                                         !! NULL.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                  optional, intent(in) :: u_accel_bt !< The zonal acceleration from terms in the
                                                     !! barotropic solver,in m s-2.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                  optional, intent(in) :: v_accel_bt !< The meridional acceleration from terms in
                                                     !! the barotropic solver,in m s-2.
  logical, optional,        intent(in) :: symmetric

!   This subroutine writes out chksums for the model's accelerations.
! Arguments: mesg - A message that appears on the chksum lines.
!  (in)      CAu - Zonal acceleration due to Coriolis and momentum
!                  advection terms, in m s-2.
!  (in)      CAv - Meridional acceleration due to Coriolis and
!                  momentum advection terms, in m s-2.
!  (in)      PFu - Zonal acceleration due to pressure gradients
!                  (equal to -dM/dx) in m s-2.
!  (in)      PFv - Meridional acceleration due to pressure
!                  gradients (equal to -dM/dy) in m s-2.
!  (in)      diffu - Zonal acceleration due to convergence of the
!                    along-isopycnal stress tensor, in m s-2.
!  (in)      diffv - Meridional acceleration due to convergence of
!                    the along-isopycnal stress tensor, in m s-2.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      pbce - the baroclinic pressure anomaly in each layer
!                   due to free surface height anomalies, in m2 s-2 H-1.
!                   pbce points to a space with nz layers or NULL.
!  (in)      u_accel_bt - The zonal acceleration from terms in the barotropic
!                         solver, in m s-2.
!  (in)      v_accel_bt - The meridional acceleration from terms in the
!                         barotropic solver, in m s-2.
  integer :: is, ie, js, je, nz
  logical :: sym

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  sym=.false.; if (present(symmetric)) sym=symmetric

  ! Note that for the chksum calls to be useful for reproducing across PE
  ! counts, there must be no redundant points, so all variables use is..ie
  ! and js...je as their extent.
  call uvchksum(mesg//" CA[uv]", cau, cav, g%HI, haloshift=0, symmetric=sym)
  call uvchksum(mesg//" PF[uv]", pfu, pfv, g%HI, haloshift=0, symmetric=sym)
  call uvchksum(mesg//" diffu", diffu, diffv, g%HI,haloshift=0, symmetric=sym)
  if (present(pbce)) &
    call hchksum(pbce, mesg//" pbce",g%HI,haloshift=0, scale=gv%m_to_H)
  if (present(u_accel_bt) .and. present(v_accel_bt)) &
    call uvchksum(mesg//" [uv]_accel_bt", u_accel_bt, v_accel_bt, g%HI,haloshift=0, symmetric=sym)
end subroutine mom_accel_chksum

! =============================================================================

subroutine mom_state_stats(mesg, u, v, h, Temp, Salt, G, allowChange, permitDiminishing)
  type(ocean_grid_type),   intent(in) :: G    !< The ocean's grid structure.
  character(len=*),        intent(in) :: mesg !< A message that appears on the chksum lines.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                           intent(in) :: u    !< The zonal velocity, in m s-1.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                           intent(in) :: v    !< The meridional velocity, in m s-1.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  &
                           intent(in) :: h    !< Layer thicknesses, in H (usually m or kg m-2).
  real, pointer, dimension(:,:,:),           &
                           intent(in) :: Temp !< Temperature in degree C.
  real, pointer, dimension(:,:,:),           &
                           intent(in) :: Salt !< Salinity, in ppt.

  logical, optional,       intent(in) :: allowChange !< do not flag an error
                                                                       !! if the statistics change.
  logical, optional,                         &
                           intent(in) :: permitDiminishing !< do not flag error
                                                           !!if the extrema are diminishing.
!   This subroutine monitors statistics for the model's state variables.
! Arguments: mesg - A message that appears on the chksum lines.
!  (in) u - Zonal velocity, in m s-1.
!  (in) v - Meridional velocity, in m s-1.
!  (in) h - Layer thickness, in m.
!  (in) T - Temperature, in degree C.
!  (in) S - Salinity, in ppt.
!  (in) G - The ocean's grid structure.
!  (in) allowChange - do not flag an error if the statistics change
!  (in) permitDiminishing - do not flag an error if the extrema are diminishing
  integer :: is, ie, js, je, nz, i, j, k
  real :: Vol, dV, Area, h_minimum
  type(stats) :: T, S, delT, delS
  type(stats), save :: oldT, oldS     ! NOTE: save data is not normally allowed but
  logical, save :: firstCall = .true. ! we use it for debugging purposes here on the
  logical :: do_TS
  real, save :: oldVol                ! assumption we will not turn this on with threads
  character(len=80) :: lMsg
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke

  do_ts = associated(temp) .and. associated(salt)

  ! First collect local stats
  area = 0. ; vol = 0.
  do j = js, je ; do i = is, ie
    area = area + g%areaT(i,j)
  enddo ; enddo
  t%minimum = 1.e34 ; t%maximum = -1.e34 ; t%average = 0.
  s%minimum = 1.e34 ; s%maximum = -1.e34 ; s%average = 0.
  h_minimum = 1.e34
  do k = 1, nz ; do j = js, je ; do i = is, ie
    if (g%mask2dT(i,j)>0.) then
      dv = g%areaT(i,j)*h(i,j,k) ; vol = vol + dv
      if (do_ts .and. h(i,j,k)>0.) then
        t%minimum = min( t%minimum, temp(i,j,k) ) ; t%maximum = max( t%maximum, temp(i,j,k) )
        t%average = t%average + dv*temp(i,j,k)
        s%minimum = min( s%minimum, salt(i,j,k) ) ; s%maximum = max( s%maximum, salt(i,j,k) )
        s%average = s%average + dv*salt(i,j,k)
      endif
      if (h_minimum > h(i,j,k)) h_minimum = h(i,j,k)
    endif
  enddo ; enddo ; enddo
  call sum_across_pes( area ) ; call sum_across_pes( vol )
  if (do_ts) then
    call min_across_pes( t%minimum ) ; call max_across_pes( t%maximum ) ; call sum_across_pes( t%average )
    call min_across_pes( s%minimum ) ; call max_across_pes( s%maximum ) ; call sum_across_pes( s%average )
    t%average = t%average / vol ; s%average = s%average / vol
  endif
  if (is_root_pe()) then
    if (.not.firstcall) then
      dv = vol - oldvol
      delt%minimum = t%minimum - oldt%minimum ; delt%maximum = t%maximum - oldt%maximum
      delt%average = t%average - oldt%average
      dels%minimum = s%minimum - olds%minimum ; dels%maximum = s%maximum - olds%maximum
      dels%average = s%average - olds%average
      write(lmsg(1:80),'(2(a,es12.4))') 'Mean thickness =',vol/area,' frac. delta=',dv/vol
      call mom_mesg(lmsg//trim(mesg))
      if (do_ts) then
        write(lmsg(1:80),'(a,3es12.4)') 'Temp min/mean/max =',t%minimum,t%average,t%maximum
        call mom_mesg(lmsg//trim(mesg))
        write(lmsg(1:80),'(a,3es12.4)') 'delT min/mean/max =',delt%minimum,delt%average,delt%maximum
        call mom_mesg(lmsg//trim(mesg))
        write(lmsg(1:80),'(a,3es12.4)') 'Salt min/mean/max =',s%minimum,s%average,s%maximum
        call mom_mesg(lmsg//trim(mesg))
        write(lmsg(1:80),'(a,3es12.4)') 'delS min/mean/max =',dels%minimum,dels%average,dels%maximum
        call mom_mesg(lmsg//trim(mesg))
      endif
    else
      write(lmsg(1:80),'(a,es12.4)') 'Mean thickness =',vol/area
      call mom_mesg(lmsg//trim(mesg))
      if (do_ts) then
        write(lmsg(1:80),'(a,3es12.4)') 'Temp min/mean/max =',t%minimum,t%average,t%maximum
        call mom_mesg(lmsg//trim(mesg))
        write(lmsg(1:80),'(a,3es12.4)') 'Salt min/mean/max =',s%minimum,s%average,s%maximum
        call mom_mesg(lmsg//trim(mesg))
      endif
    endif
  endif
  firstcall = .false. ; oldvol = vol
  oldt%minimum = t%minimum ; oldt%maximum = t%maximum ; oldt%average = t%average
  olds%minimum = s%minimum ; olds%maximum = s%maximum ; olds%average = s%average

  if (do_ts .and. t%minimum<-5.0) then
    do j = js, je ; do i = is, ie
      if (minval(temp(i,j,:)) == t%minimum) then
        write(0,'(a,2f12.5)') 'x,y=',g%geoLonT(i,j),g%geoLatT(i,j)
        write(0,'(a3,3a12)') 'k','h','Temp','Salt'
        do k = 1, nz
          write(0,'(i3,3es12.4)') k,h(i,j,k),temp(i,j,k),salt(i,j,k)
        enddo
        stop 'Extremum detected'
      endif
    enddo ; enddo
  endif

  if (h_minimum<0.0) then
    do j = js, je ; do i = is, ie
      if (minval(h(i,j,:)) == h_minimum) then
        write(0,'(a,2f12.5)') 'x,y=',g%geoLonT(i,j),g%geoLatT(i,j)
        write(0,'(a3,3a12)') 'k','h','Temp','Salt'
        do k = 1, nz
          write(0,'(i3,3es12.4)') k,h(i,j,k),temp(i,j,k),salt(i,j,k)
        enddo
        stop 'Negative thickness detected'
      endif
    enddo ; enddo
  endif

end subroutine mom_state_stats

end module mom_checksum_packages