MOM_PointAccel.F90

!> Debug accelerations at a given point
!!
!!    The two subroutines in this file write out all of the terms
!! in the u- or v-momentum balance at a given point.  Usually
!! these subroutines are called after the velocities exceed some
!! threshold, in order to determine which term is culpable.
!! often this is done for debugging purposes.
module mom_pointaccel

! This file is part of MOM6. See LICENSE.md for the license.

use mom_diag_mediator, only : diag_ctrl
use mom_domains, only : pe_here
use mom_error_handler, only : mom_error, note
use mom_file_parser, only : get_param, log_param, log_version, param_file_type
use mom_get_input, only : directories
use mom_grid, only : ocean_grid_type
use mom_io, only : open_file
use mom_io, only : append_file, ascii_file, multiple, single_file
use mom_time_manager, only : time_type, get_time, get_date, set_date, operator(-)
use mom_unit_scaling, only : unit_scale_type
use mom_variables, only : ocean_internal_state, accel_diag_ptrs, cont_diag_ptrs
use mom_verticalgrid, only : verticalgrid_type

implicit none ; private

#include <MOM_memory.h>

public write_u_accel, write_v_accel, pointaccel_init

!> The control structure for the MOM_PointAccel module
type, public :: pointaccel_cs ; private
  character(len=200) :: u_trunc_file !< The complete path to the file in which a column's worth of
                                     !! u-accelerations are written if u-velocity truncations occur.
  character(len=200) :: v_trunc_file !< The complete path to the file in which a column's worth of
                                     !! v-accelerations are written if v-velocity truncations occur.
  integer :: u_file         !< The unit number for an opened u-truncation files, or -1 if it has not yet been opened.
  integer :: v_file         !< The unit number for an opened v-truncation files, or -1 if it has not yet been opened.
  integer :: cols_written   !< The number of columns whose output has been
                            !! written by this PE during the current run.
  integer :: max_writes     !< The maximum number of times any PE can write out
                            !! a column's worth of accelerations during a run.
  type(time_type), pointer :: time => null() !< A pointer to the ocean model's clock.
  type(diag_ctrl), pointer :: diag => null() !< A structure that is used to
                                   !! regulate the timing of diagnostic output.
! The following are pointers to many of the state variables and accelerations
! that are used to step the physical model forward.  They all use the same
! names as the variables they point to in MOM.F90
  real, pointer, dimension(:,:,:) :: &
    u_av => null(), &       !< Time average u-velocity [L T-1 ~> m s-1].
    v_av => null(), &       !< Time average velocity [L T-1 ~> m s-1].
    u_prev => null(), &     !< Previous u-velocity [L T-1 ~> m s-1].
    v_prev => null(), &     !< Previous v-velocity [L T-1 ~> m s-1].
    t => null(), &          !< Temperature [degC].
    s => null(), &          !< Salinity [ppt].
    u_accel_bt => null(), & !< Barotropic u-acclerations [L T-2 ~> m s-2]
    v_accel_bt => null()    !< Barotropic v-acclerations [L T-2 ~> m s-2]
  real, pointer, dimension(:,:,:) :: pbce => null() !< pbce times eta gives the baroclinic
                            !! pressure anomaly in each layer due to free surface height anomalies
                            !! [m2 s-2 H-1 ~> m s-2 or m4 kg-1 s-2].
end type pointaccel_cs

contains

!> This subroutine writes to an output file all of the accelerations
!! that have been applied to a column of zonal velocities over the
!! previous timestep.  This subroutine is called from vertvisc.
subroutine write_u_accel(I, j, um, hin, ADp, CDp, dt_in_T, G, GV, US, CS, vel_rpt, str, a, hv)
  integer,                     intent(in) :: i   !< The zonal index of the column to be documented.
  integer,                     intent(in) :: j   !< The meridional index of the column to be documented.
  type(ocean_grid_type),       intent(in) :: g   !< The ocean's grid structure.
  type(verticalgrid_type),     intent(in) :: gv  !< The ocean's vertical grid structure.
  type(unit_scale_type),       intent(in) :: us  !< A dimensional unit scaling type
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                               intent(in) :: um  !< The new zonal velocity [L T-1 ~> m s-1].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  &
                               intent(in) :: hin !< The layer thickness [H ~> m or kg m-2].
  type(accel_diag_ptrs),       intent(in) :: adp !< A structure pointing to the various
                                                 !! accelerations in the momentum equations.
  type(cont_diag_ptrs),        intent(in) :: cdp !<  A structure with pointers to various terms
                                                 !! in the continuity equations.
  real,                        intent(in) :: dt_in_t  !< The ocean dynamics time step [T ~> s].
  type(pointaccel_cs),         pointer    :: cs  !< The control structure returned by a previous
                                                 !! call to PointAccel_init.
  real,                        intent(in) :: vel_rpt !< The velocity magnitude that triggers a report [L T-1 ~> m s-1].
  real, optional,              intent(in) :: str !< The surface wind stress integrated over a time
                                                 !! step divided by the Boussinesq density [m2 s-1].
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                     optional, intent(in) :: a   !< The layer coupling coefficients from vertvisc [Z s-1 ~> m s-1].
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                     optional, intent(in) :: hv  !< The layer thicknesses at velocity grid points,
                                                 !! from vertvisc [H ~> m or kg m-2].
  ! Local variables
  real    :: f_eff, cfl
  real    :: angstrom
  real    :: truncvel, du
  real    :: dt  ! The time step [s]
  real    :: inorm(szk_(g))
  real    :: e(szk_(g)+1)
  real    :: h_scale, uh_scale
  integer :: yr, mo, day, hr, minute, sec, yearday
  integer :: k, ks, ke
  integer :: nz
  logical :: do_k(szk_(g)+1)
  logical :: prev_avail
  integer :: file

  angstrom = gv%Angstrom_H + gv%H_subroundoff
  dt = us%T_to_s*dt_in_t
  h_scale = gv%H_to_m ; uh_scale = gv%H_to_m*us%L_T_to_m_s

!  if (.not.associated(CS)) return
  nz = g%ke
  if (cs%cols_written < cs%max_writes) then
    cs%cols_written = cs%cols_written + 1

    ks = 1 ; ke = nz
    do_k(:) = .false.

  ! Open up the file for output if this is the first call.
    if (cs%u_file < 0) then
      if (len_trim(cs%u_trunc_file) < 1) return
      call open_file(cs%u_file, trim(cs%u_trunc_file), action=append_file, &
                     form=ascii_file, threading=multiple, fileset=single_file)
      if (cs%u_file < 0) then
        call mom_error(note, 'Unable to open file '//trim(cs%u_trunc_file)//'.')
        return
      endif
    endif
    file = cs%u_file

    prev_avail = (associated(cs%u_prev) .and. associated(cs%v_prev))

  ! Determine which layers to write out accelerations for.
    do k=1,nz
      if (((max(cs%u_av(i,j,k),um(i,j,k)) >= vel_rpt) .or. &
           (min(cs%u_av(i,j,k),um(i,j,k)) <= -vel_rpt)) .and. &
          ((hin(i,j,k) + hin(i+1,j,k)) > 3.0*angstrom)) exit
    enddo
    ks = k
    do k=nz,1,-1
      if (((max(cs%u_av(i,j,k), um(i,j,k)) >= vel_rpt) .or. &
           (min(cs%u_av(i,j,k), um(i,j,k)) <= -vel_rpt)) .and. &
          ((hin(i,j,k) + hin(i+1,j,k)) > 3.0*angstrom)) exit
    enddo
    ke = k
    if (ke < ks) then
      ks = 1; ke = nz; write(file,'("U: Unable to set ks & ke.")')
    endif

    call get_date(cs%Time, yr, mo, day, hr, minute, sec)
    call get_time((cs%Time - set_date(yr, 1, 1, 0, 0, 0)), sec, yearday)
    write (file,'(/,"--------------------------")')
    write (file,'(/,"Time ",i5,i4,F6.2," U-velocity violation at ",I4,": ",2(I3), &
        & " (",F7.2," E "F7.2," N) Layers ",I3," to ",I3,". dt = ",1PG10.4)') &
        yr, yearday, (real(sec)/3600.0), pe_here(), i, j, &
        g%geoLonCu(i,j), g%geoLatCu(i,j), ks, ke, dt

    if (ks <= gv%nk_rho_varies) ks = 1
    do k=ks,ke
      if ((hin(i,j,k) + hin(i+1,j,k)) > 3.0*angstrom) do_k(k) = .true.
    enddo

    write(file,'(/,"Layers:",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(I10," ",$)') (k); enddo
    write(file,'(/,"u(m):  ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (us%L_T_to_m_s*um(i,j,k)); enddo
    if (prev_avail) then
      write(file,'(/,"u(mp): ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (us%L_T_to_m_s*cs%u_prev(i,j,k)); enddo
    endif
    write(file,'(/,"u(3):  ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (us%L_T_to_m_s*cs%u_av(i,j,k)); enddo

    write(file,'(/,"CFL u: ",$)')
    do k=ks,ke ; if (do_k(k)) then
      cfl = abs(um(i,j,k)) * us%s_to_T*dt * g%dy_Cu(i,j)
      if (um(i,j,k) < 0.0) then ; cfl = cfl * g%IareaT(i+1,j)
      else ; cfl = cfl * g%IareaT(i,j) ; endif
      write(file,'(ES10.3," ",$)') cfl
    endif ; enddo
    write(file,'(/,"CFL0 u:",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                    abs(um(i,j,k)) * us%s_to_T*dt * g%IdxCu(i,j) ; enddo

    if (prev_avail) then
      write(file,'(/,"du:    ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                      (us%L_T_to_m_s*(um(i,j,k)-cs%u_prev(i,j,k))); enddo
    endif
    write(file,'(/,"CAu:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (dt*us%L_T2_to_m_s2*adp%CAu(i,j,k)); enddo
    write(file,'(/,"PFu:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (dt*us%L_T2_to_m_s2*adp%PFu(i,j,k)); enddo
    write(file,'(/,"diffu: ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (dt*us%L_T2_to_m_s2*adp%diffu(i,j,k)); enddo

    if (associated(adp%gradKEu)) then
      write(file,'(/,"KEu:   ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                      (dt*us%L_T2_to_m_s2*adp%gradKEu(i,j,k)); enddo
    endif
    if (associated(adp%rv_x_v)) then
      write(file,'(/,"Coru:  ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
          dt*us%L_T2_to_m_s2*(adp%CAu(i,j,k)-adp%rv_x_v(i,j,k)); enddo
    endif
    if (associated(adp%du_dt_visc)) then
      write(file,'(/,"ubv:   ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
          us%L_T_to_m_s*(um(i,j,k) - us%s_to_T*dt*adp%du_dt_visc(i,j,k)); enddo
      write(file,'(/,"duv:   ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                      (dt*us%L_T2_to_m_s2*adp%du_dt_visc(i,j,k)); enddo
    endif
    if (associated(adp%du_other)) then
      write(file,'(/,"du_other: ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                      (us%L_T_to_m_s*adp%du_other(i,j,k)); enddo
    endif
    if (present(a)) then
      write(file,'(/,"a:     ",$)')
      do k=ks,ke+1 ; if (do_k(k)) write(file,'(ES10.3," ",$)') a(i,j,k)*us%Z_to_m*dt; enddo
    endif
    if (present(hv)) then
      write(file,'(/,"hvel:  ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') hv(i,j,k); enddo
    endif
    write(file,'(/,"Stress:  ",ES10.3)') str

    if (associated(cs%u_accel_bt)) then
      write(file,'("dubt:  ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                      (dt*us%L_T2_to_m_s2*cs%u_accel_bt(i,j,k)) ; enddo
      write(file,'(/)')
    endif

    write(file,'(/,"h--:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (h_scale*hin(i,j-1,k)); enddo
    write(file,'(/,"h+-:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (h_scale*hin(i+1,j-1,k)); enddo
    write(file,'(/,"h-0:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (h_scale*hin(i,j,k)); enddo
    write(file,'(/,"h+0:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (h_scale*hin(i+1,j,k)); enddo
    write(file,'(/,"h-+:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (h_scale*hin(i,j+1,k)); enddo
    write(file,'(/,"h++:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (h_scale*hin(i+1,j+1,k)); enddo


    e(nz+1) = -us%Z_to_m*g%bathyT(i,j)
    do k=nz,1,-1 ; e(k) = e(k+1) + h_scale*hin(i,j,k) ; enddo
    write(file,'(/,"e-:    ",$)')
    write(file,'(ES10.3," ",$)') e(ks)
    do k=ks+1,ke+1 ; if (do_k(k-1)) write(file,'(ES10.3," ",$)') e(k); enddo

    e(nz+1) = -us%Z_to_m*g%bathyT(i+1,j)
    do k=nz,1,-1 ; e(k) = e(k+1) + h_scale*hin(i+1,j,k) ; enddo
    write(file,'(/,"e+:    ",$)')
    write(file,'(ES10.3," ",$)') e(ks)
    do k=ks+1,ke+1 ; if (do_k(k-1)) write(file,'(ES10.3," ",$)') e(k) ; enddo
    if (associated(cs%T)) then
      write(file,'(/,"T-:    ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') cs%T(i,j,k); enddo
      write(file,'(/,"T+:    ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') cs%T(i+1,j,k); enddo
    endif
    if (associated(cs%S)) then
      write(file,'(/,"S-:    ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') cs%S(i,j,k); enddo
      write(file,'(/,"S+:    ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') cs%S(i+1,j,k); enddo
    endif

    if (prev_avail) then
      write(file,'(/,"v--:   ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (cs%v_prev(i,j-1,k)); enddo
      write(file,'(/,"v-+:   ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (cs%v_prev(i,j,k)); enddo
      write(file,'(/,"v+-:   ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (cs%v_prev(i+1,j-1,k)); enddo
      write(file,'(/,"v++:   ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (cs%v_prev(i+1,j,k)); enddo
    endif

    write(file,'(/,"vh--:  ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                    (uh_scale*cdp%vh(i,j-1,k)*g%IdxCv(i,j-1)); enddo
    write(file,'(/," vhC--:",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                        (0.5*us%L_T_to_m_s*cs%v_av(i,j-1,k)*h_scale*(hin(i,j-1,k) + hin(i,j,k))); enddo
    if (prev_avail) then
      write(file,'(/," vhCp--:",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                          (0.5*cs%v_prev(i,j-1,k)*h_scale*(hin(i,j-1,k) + hin(i,j,k))); enddo
    endif

    write(file,'(/,"vh-+:  ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                    (uh_scale*cdp%vh(i,j,k)*g%IdxCv(i,j)); enddo
    write(file,'(/," vhC-+:",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                        (0.5*us%L_T_to_m_s*cs%v_av(i,j,k)*h_scale*(hin(i,j,k) + hin(i,j+1,k))); enddo
    if (prev_avail) then
      write(file,'(/," vhCp-+:",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                          (0.5*cs%v_prev(i,j,k)*h_scale*(hin(i,j,k) + hin(i,j+1,k))); enddo
    endif

    write(file,'(/,"vh+-:  ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                      (uh_scale*cdp%vh(i+1,j-1,k)*g%IdxCv(i+1,j-1)); enddo
    write(file,'(/," vhC+-:",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                    (0.5*us%L_T_to_m_s*cs%v_av(i+1,j-1,k)*h_scale*(hin(i+1,j-1,k) + hin(i+1,j,k))); enddo
    if (prev_avail) then
      write(file,'(/," vhCp+-:",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                      (0.5*cs%v_prev(i+1,j-1,k)*h_scale*(hin(i+1,j-1,k) + hin(i+1,j,k))); enddo
    endif

    write(file,'(/,"vh++:  ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                          (uh_scale*cdp%vh(i+1,j,k)*g%IdxCv(i+1,j)); enddo
    write(file,'(/," vhC++:",$)')
         do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                     (0.5*us%L_T_to_m_s*cs%v_av(i+1,j,k)*h_scale*(hin(i+1,j,k) + hin(i+1,j+1,k))); enddo
    if (prev_avail) then
      write(file,'(/," vhCp++:",$)')
           do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                       (0.5*us%L_T_to_m_s*cs%v_av(i+1,j,k)*h_scale*(hin(i+1,j,k) + hin(i+1,j+1,k))); enddo
    endif

    write(file,'(/,"D:     ",2(ES10.3))') us%Z_to_m*g%bathyT(i,j),us%Z_to_m*g%bathyT(i+1,j)

  !  From here on, the normalized accelerations are written.
    if (prev_avail) then
      do k=ks,ke
        du = us%L_T_to_m_s*(um(i,j,k) - cs%u_prev(i,j,k))
        if (abs(du) < 1.0e-6) du = 1.0e-6
        inorm(k) = 1.0 / du
      enddo

      write(file,'(2/,"Norm:  ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') (1.0/inorm(k)); enddo

      write(file,'(/,"du:    ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') &
                        (us%L_T_to_m_s*(um(i,j,k)-cs%u_prev(i,j,k))*inorm(k)); enddo

      write(file,'(/,"CAu:   ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') &
                                      (dt*us%L_T2_to_m_s2*adp%CAu(i,j,k)*inorm(k)); enddo

      write(file,'(/,"PFu:   ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') &
                                      (dt*us%L_T2_to_m_s2*adp%PFu(i,j,k)*inorm(k)); enddo

      write(file,'(/,"diffu: ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') &
                                      (dt*us%L_T2_to_m_s2*adp%diffu(i,j,k)*inorm(k)); enddo

      if (associated(adp%gradKEu)) then
        write(file,'(/,"KEu:   ",$)')
        do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') &
                                        (dt*us%L_T2_to_m_s2*adp%gradKEu(i,j,k)*inorm(k)); enddo
      endif
      if (associated(adp%rv_x_v)) then
        write(file,'(/,"Coru:  ",$)')
        do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') &
            dt*us%L_T2_to_m_s2*(adp%CAu(i,j,k)-adp%rv_x_v(i,j,k))*inorm(k); enddo
      endif
      if (associated(adp%du_dt_visc)) then
        write(file,'(/,"duv:   ",$)')
        do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') &
            (dt*us%L_T2_to_m_s2*adp%du_dt_visc(i,j,k))*inorm(k); enddo
      endif
      if (associated(adp%du_other)) then
        write(file,'(/,"du_other: ",$)')
        do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') &
            (us%L_T_to_m_s*adp%du_other(i,j,k))*inorm(k); enddo
      endif
      if (associated(cs%u_accel_bt)) then
        write(file,'(/,"dubt:  ",$)')
        do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') &
                                        (dt*us%L_T2_to_m_s2*cs%u_accel_bt(i,j,k)*inorm(k)) ; enddo
      endif
    endif

    write(file,'(2/)')

    call flush(file)
  endif

end subroutine write_u_accel

!> This subroutine writes to an output file all of the accelerations
!! that have been applied to a column of meridional velocities over
!! the previous timestep.  This subroutine is called from vertvisc.
subroutine write_v_accel(i, J, vm, hin, ADp, CDp, dt_in_T, G, GV, US, CS, vel_rpt, str, a, hv)
  integer,                     intent(in) :: i   !< The zonal index of the column to be documented.
  integer,                     intent(in) :: j   !< The meridional index of the column to be documented.
  type(ocean_grid_type),       intent(in) :: g   !< The ocean's grid structure.
  type(verticalgrid_type),     intent(in) :: gv  !< The ocean's vertical grid structure.
  type(unit_scale_type),       intent(in) :: us  !< A dimensional unit scaling type
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                               intent(in) :: vm  !< The new meridional velocity [L T-1 ~> m s-1].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  &
                               intent(in) :: hin !< The layer thickness [H ~> m or kg m-2].
  type(accel_diag_ptrs),       intent(in) :: adp !< A structure pointing to the various
                                                 !! accelerations in the momentum equations.
  type(cont_diag_ptrs),        intent(in) :: cdp !< A structure with pointers to various terms in
                                                 !! the continuity equations.
  real,                        intent(in) :: dt_in_t  !< The ocean dynamics time step [T ~> s].
  type(pointaccel_cs),         pointer    :: cs  !< The control structure returned by a previous
                                                 !! call to PointAccel_init.
  real,                        intent(in) :: vel_rpt !< The velocity magnitude that triggers a report [L T-1 ~> m s-1].
  real, optional,              intent(in) :: str !< The surface wind stress integrated over a time
                                                 !! step divided by the Boussinesq density [m2 s-1].
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                     optional, intent(in) :: a   !< The layer coupling coefficients from vertvisc [Z s-1 ~> m s-1].
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                     optional, intent(in) :: hv  !< The layer thicknesses at velocity grid points,
                                                 !! from vertvisc [H ~> m or kg m-2].
  ! Local variables
  real    :: f_eff, cfl
  real    :: angstrom
  real    :: truncvel, dv
  real    :: dt  ! The time step [s]
  real    :: inorm(szk_(g))
  real    :: e(szk_(g)+1)
  real    :: h_scale, uh_scale
  integer :: yr, mo, day, hr, minute, sec, yearday
  integer :: k, ks, ke
  integer :: nz
  logical :: do_k(szk_(g)+1)
  logical :: prev_avail
  integer :: file

  angstrom = gv%Angstrom_H + gv%H_subroundoff
  dt = us%T_to_s*dt_in_t
  h_scale = gv%H_to_m ; uh_scale = gv%H_to_m*us%L_T_to_m_s

!  if (.not.associated(CS)) return
  nz = g%ke
  if (cs%cols_written < cs%max_writes) then
    cs%cols_written = cs%cols_written + 1

    ks = 1 ; ke = nz
    do_k(:) = .false.

  ! Open up the file for output if this is the first call.
    if (cs%v_file < 0) then
      if (len_trim(cs%v_trunc_file) < 1) return
      call open_file(cs%v_file, trim(cs%v_trunc_file), action=append_file, &
                     form=ascii_file, threading=multiple, fileset=single_file)
      if (cs%v_file < 0) then
        call mom_error(note, 'Unable to open file '//trim(cs%v_trunc_file)//'.')
        return
      endif
    endif
    file = cs%v_file

    prev_avail = (associated(cs%u_prev) .and. associated(cs%v_prev))

    do k=1,nz
      if (((max(cs%v_av(i,j,k), vm(i,j,k)) >= vel_rpt) .or. &
           (min(cs%v_av(i,j,k), vm(i,j,k)) <= -vel_rpt)) .and. &
          ((hin(i,j,k) + hin(i,j+1,k)) > 3.0*angstrom)) exit
    enddo
    ks = k
    do k=nz,1,-1
      if (((max(cs%v_av(i,j,k), vm(i,j,k)) >= vel_rpt) .or. &
           (min(cs%v_av(i,j,k), vm(i,j,k)) <= -vel_rpt)) .and. &
          ((hin(i,j,k) + hin(i,j+1,k)) > 3.0*angstrom)) exit
    enddo
    ke = k
    if (ke < ks) then
      ks = 1; ke = nz; write(file,'("V: Unable to set ks & ke.")')
    endif

    call get_date(cs%Time, yr, mo, day, hr, minute, sec)
    call get_time((cs%Time - set_date(yr, 1, 1, 0, 0, 0)), sec, yearday)
    write (file,'(/,"--------------------------")')
    write (file,'(/,"Time ",i5,i4,F6.2," V-velocity violation at ",I4,": ",2(I3), &
        & " (",F7.2," E ",F7.2," N) Layers ",I3," to ",I3,". dt = ",1PG10.4)') &
        yr, yearday, (real(sec)/3600.0), pe_here(), i, j, &
        g%geoLonCv(i,j), g%geoLatCv(i,j), ks, ke, dt

    if (ks <= gv%nk_rho_varies) ks = 1
    do k=ks,ke
      if ((hin(i,j,k) + hin(i,j+1,k)) > 3.0*angstrom) do_k(k) = .true.
    enddo

    write(file,'(/,"Layers:",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(I10," ",$)') (k); enddo
    write(file,'(/,"v(m):  ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (us%L_T_to_m_s*vm(i,j,k)); enddo

    if (prev_avail) then
      write(file,'(/,"v(mp): ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (us%L_T_to_m_s*cs%v_prev(i,j,k)); enddo
    endif

    write(file,'(/,"v(3):  ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (us%L_T_to_m_s*cs%v_av(i,j,k)); enddo
    write(file,'(/,"CFL v: ",$)')
    do k=ks,ke ; if (do_k(k)) then
      cfl = abs(vm(i,j,k)) * us%s_to_T*dt * g%dx_Cv(i,j)
      if (vm(i,j,k) < 0.0) then ; cfl = cfl * g%IareaT(i,j+1)
      else ; cfl = cfl * g%IareaT(i,j) ; endif
      write(file,'(ES10.3," ",$)') cfl
    endif ; enddo
    write(file,'(/,"CFL0 v:",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                    abs(vm(i,j,k)) * us%s_to_T*dt * g%IdyCv(i,j) ; enddo

    if (prev_avail) then
      write(file,'(/,"dv:    ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                      (us%L_T_to_m_s*(vm(i,j,k)-cs%v_prev(i,j,k))); enddo
    endif

    write(file,'(/,"CAv:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (dt*us%L_T2_to_m_s2*adp%CAv(i,j,k)); enddo

    write(file,'(/,"PFv:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (dt*us%L_T2_to_m_s2*adp%PFv(i,j,k)); enddo

    write(file,'(/,"diffv: ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (dt*us%L_T2_to_m_s2*adp%diffv(i,j,k)); enddo

    if (associated(adp%gradKEv)) then
      write(file,'(/,"KEv:   ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                      (dt*us%L_T2_to_m_s2*adp%gradKEv(i,j,k)); enddo
    endif
    if (associated(adp%rv_x_u)) then
      write(file,'(/,"Corv:  ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                 dt*us%L_T2_to_m_s2*(adp%CAv(i,j,k)-adp%rv_x_u(i,j,k)); enddo
    endif
    if (associated(adp%dv_dt_visc)) then
      write(file,'(/,"vbv:   ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
          us%L_T_to_m_s*(vm(i,j,k) - us%s_to_T*dt*adp%dv_dt_visc(i,j,k)); enddo

      write(file,'(/,"dvv:   ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                      (dt*us%L_T2_to_m_s2*adp%dv_dt_visc(i,j,k)); enddo
    endif
    if (associated(adp%dv_other)) then
      write(file,'(/,"dv_other: ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                      (us%L_T_to_m_s*adp%dv_other(i,j,k)); enddo
    endif
    if (present(a)) then
      write(file,'(/,"a:     ",$)')
      do k=ks,ke+1 ; if (do_k(k)) write(file,'(ES10.3," ",$)') a(i,j,k)*us%Z_to_m*dt; enddo
    endif
    if (present(hv)) then
      write(file,'(/,"hvel:  ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') hv(i,j,k); enddo
    endif
    write(file,'(/,"Stress:  ",ES10.3)') str

    if (associated(cs%v_accel_bt)) then
      write(file,'("dvbt:  ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                      (dt*us%L_T2_to_m_s2*cs%v_accel_bt(i,j,k)) ; enddo
      write(file,'(/)')
    endif

    write(file,'("h--:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') h_scale*hin(i-1,j,k); enddo
    write(file,'(/,"h0-:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') h_scale*hin(i,j,k); enddo
    write(file,'(/,"h+-:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') h_scale*hin(i+1,j,k); enddo
    write(file,'(/,"h-+:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') h_scale*hin(i-1,j+1,k); enddo
    write(file,'(/,"h0+:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') h_scale*hin(i,j+1,k); enddo
    write(file,'(/,"h++:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') h_scale*hin(i+1,j+1,k); enddo

    e(nz+1) = -us%Z_to_m*g%bathyT(i,j)
    do k=nz,1,-1 ; e(k) = e(k+1) + h_scale*hin(i,j,k); enddo
    write(file,'(/,"e-:    ",$)')
    write(file,'(ES10.3," ",$)') e(ks)
    do k=ks+1,ke+1 ; if (do_k(k-1)) write(file,'(ES10.3," ",$)') e(k); enddo

    e(nz+1) = -us%Z_to_m*g%bathyT(i,j+1)
    do k=nz,1,-1 ; e(k) = e(k+1) + h_scale*hin(i,j+1,k) ; enddo
    write(file,'(/,"e+:    ",$)')
    write(file,'(ES10.3," ",$)') e(ks)
    do k=ks+1,ke+1 ; if (do_k(k-1)) write(file,'(ES10.3," ",$)') e(k); enddo
    if (associated(cs%T)) then
      write(file,'(/,"T-:    ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') cs%T(i,j,k); enddo
      write(file,'(/,"T+:    ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') cs%T(i,j+1,k); enddo
    endif
    if (associated(cs%S)) then
      write(file,'(/,"S-:    ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') cs%S(i,j,k); enddo
      write(file,'(/,"S+:    ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') cs%S(i,j+1,k); enddo
    endif

    if (prev_avail) then
      write(file,'(/,"u--:   ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') cs%u_prev(i-1,j,k); enddo
      write(file,'(/,"u-+:   ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') cs%u_prev(i-1,j+1,k); enddo
      write(file,'(/,"u+-:   ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') cs%u_prev(i,j,k); enddo
      write(file,'(/,"u++:   ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') cs%u_prev(i,j+1,k); enddo
    endif

    write(file,'(/,"uh--:  ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                    (uh_scale*cdp%uh(i-1,j,k)*g%IdyCu(i-1,j)); enddo
    write(file,'(/," uhC--: ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
            (us%L_T_to_m_s*cs%u_av(i-1,j,k) * h_scale*0.5*(hin(i-1,j,k) + hin(i,j,k))); enddo
    if (prev_avail) then
      write(file,'(/," uhCp--:",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
            (cs%u_prev(i-1,j,k) * h_scale*0.5*(hin(i-1,j,k) + hin(i,j,k))); enddo
    endif

    write(file,'(/,"uh-+:  ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                    (uh_scale*cdp%uh(i-1,j+1,k)*g%IdyCu(i-1,j+1)); enddo
    write(file,'(/," uhC-+: ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
            (us%L_T_to_m_s*cs%u_av(i-1,j+1,k) * h_scale*0.5*(hin(i-1,j+1,k) + hin(i,j+1,k))); enddo
    if (prev_avail) then
      write(file,'(/," uhCp-+:",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
            (cs%u_prev(i-1,j+1,k) * h_scale*0.5*(hin(i-1,j+1,k) + hin(i,j+1,k))); enddo
    endif

    write(file,'(/,"uh+-:  ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                    (uh_scale*cdp%uh(i,j,k)*g%IdyCu(i,j)); enddo
    write(file,'(/," uhC+-: ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
            (us%L_T_to_m_s*cs%u_av(i,j,k) * h_scale*0.5*(hin(i,j,k) + hin(i+1,j,k))); enddo
    if (prev_avail) then
      write(file,'(/," uhCp+-:",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
            (cs%u_prev(i,j,k) * h_scale*0.5*(hin(i,j,k) + hin(i+1,j,k))); enddo
    endif

    write(file,'(/,"uh++:  ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
                                    (uh_scale*cdp%uh(i,j+1,k)*g%IdyCu(i,j+1)); enddo
    write(file,'(/," uhC++: ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
            (us%L_T_to_m_s*cs%u_av(i,j+1,k) * 0.5*h_scale*(hin(i,j+1,k) + hin(i+1,j+1,k))); enddo
    if (prev_avail) then
      write(file,'(/," uhCp++:",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') &
            (cs%u_prev(i,j+1,k) * h_scale*0.5*(hin(i,j+1,k) + hin(i+1,j+1,k))); enddo
    endif

    write(file,'(/,"D:     ",2(ES10.3))') us%Z_to_m*g%bathyT(i,j),us%Z_to_m*g%bathyT(i,j+1)

  !  From here on, the normalized accelerations are written.
    if (prev_avail) then
      do k=ks,ke
        dv = us%L_T_to_m_s*(vm(i,j,k)-cs%v_prev(i,j,k))
        if (abs(dv) < 1.0e-6) dv = 1.0e-6
        inorm(k) = 1.0 / dv
      enddo

      write(file,'(2/,"Norm:  ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') (1.0/inorm(k)); enddo
      write(file,'(/,"dv:    ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') &
                                      (us%L_T_to_m_s*(vm(i,j,k)-cs%v_prev(i,j,k))*inorm(k)); enddo
      write(file,'(/,"CAv:   ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') &
                                      (dt*us%L_T2_to_m_s2*adp%CAv(i,j,k)*inorm(k)); enddo
      write(file,'(/,"PFv:   ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') &
                                      (dt*us%L_T2_to_m_s2*adp%PFv(i,j,k)*inorm(k)); enddo
      write(file,'(/,"diffv: ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') &
                                      (dt*us%L_T2_to_m_s2*adp%diffv(i,j,k)*inorm(k)); enddo

      if (associated(adp%gradKEu)) then
        write(file,'(/,"KEv:   ",$)')
        do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') &
                                      (dt*us%L_T2_to_m_s2*adp%gradKEv(i,j,k)*inorm(k)); enddo
      endif
      if (associated(adp%rv_x_u)) then
        write(file,'(/,"Corv:  ",$)')
        do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') &
            dt*us%L_T2_to_m_s2*(adp%CAv(i,j,k)-adp%rv_x_u(i,j,k))*inorm(k); enddo
      endif
      if (associated(adp%dv_dt_visc)) then
        write(file,'(/,"dvv:   ",$)')
        do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') &
            (dt*us%L_T2_to_m_s2*adp%dv_dt_visc(i,j,k)*inorm(k)); enddo
      endif
      if (associated(adp%dv_other)) then
        write(file,'(/,"dv_other: ",$)')
        do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') &
            (us%L_T_to_m_s*adp%dv_other(i,j,k)*inorm(k)); enddo
      endif
      if (associated(cs%v_accel_bt)) then
        write(file,'(/,"dvbt:  ",$)')
        do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') &
                                        (dt*us%L_T2_to_m_s2*cs%v_accel_bt(i,j,k)*inorm(k)) ; enddo
      endif
    endif

    write(file,'(2/)')

    call flush(file)
  endif

end subroutine write_v_accel

!> This subroutine initializes the parameters regulating how truncations are logged.
subroutine pointaccel_init(MIS, Time, G, param_file, diag, dirs, CS)
  type(ocean_internal_state), &
                        target, intent(in)    :: mis  !< For "MOM Internal State" a set of pointers
                                                      !! to the fields and accelerations that make
                                                      !! up the ocean's physical state.
  type(time_type),      target, intent(in)    :: time !< The current model time.
  type(ocean_grid_type),        intent(in)    :: g    !< The ocean's grid structure.
  type(param_file_type),        intent(in)    :: param_file !< A structure to parse for run-time
                                                      !! parameters.
  type(diag_ctrl),      target, intent(inout) :: diag !< A structure that is used to regulate
                                                      !! diagnostic output.
  type(directories),            intent(in)    :: dirs !< A structure containing several relevant
                                                      !! directory paths.
  type(pointaccel_cs),          pointer       :: cs   !< A pointer that is set to point to the
                                                      !! control structure for this module.
! This include declares and sets the variable "version".
#include "version_variable.h"
  character(len=40)  :: mdl = "MOM_PointAccel" ! This module's name.

  if (associated(cs)) return
  allocate(cs)

  cs%diag => diag ; cs%Time => time

  cs%T => mis%T ; cs%S => mis%S ; cs%pbce => mis%pbce
  cs%u_accel_bt => mis%u_accel_bt ; cs%v_accel_bt => mis%v_accel_bt
  cs%u_prev => mis%u_prev ; cs%v_prev => mis%v_prev
  cs%u_av => mis%u_av; if (.not.associated(mis%u_av)) cs%u_av => mis%u(:,:,:)
  cs%v_av => mis%v_av; if (.not.associated(mis%v_av)) cs%v_av => mis%v(:,:,:)

  ! Read all relevant parameters and write them to the model log.
  call log_version(param_file, mdl, version, "", debugging=.true.)
  call get_param(param_file, mdl, "U_TRUNC_FILE", cs%u_trunc_file, &
                 "The absolute path to the file where the accelerations "//&
                 "leading to zonal velocity truncations are written. \n"//&
                 "Leave this empty for efficiency if this diagnostic is "//&
                 "not needed.", default="", debuggingparam=.true.)
  call get_param(param_file, mdl, "V_TRUNC_FILE", cs%v_trunc_file, &
                 "The absolute path to the file where the accelerations "//&
                 "leading to meridional velocity truncations are written. \n"//&
                 "Leave this empty for efficiency if this diagnostic is "//&
                 "not needed.", default="", debuggingparam=.true.)
  call get_param(param_file, mdl, "MAX_TRUNC_FILE_SIZE_PER_PE", cs%max_writes, &
                 "The maximum number of columns of truncations that any PE "//&
                 "will write out during a run.", default=50, debuggingparam=.true.)

  if (len_trim(dirs%output_directory) > 0) then
    if (len_trim(cs%u_trunc_file) > 0) &
      cs%u_trunc_file = trim(dirs%output_directory)//trim(cs%u_trunc_file)
    if (len_trim(cs%v_trunc_file) > 0) &
      cs%v_trunc_file = trim(dirs%output_directory)//trim(cs%v_trunc_file)
    call log_param(param_file, mdl, "output_dir/U_TRUNC_FILE", cs%u_trunc_file, debuggingparam=.true.)
    call log_param(param_file, mdl, "output_dir/V_TRUNC_FILE", cs%v_trunc_file, debuggingparam=.true.)
  endif
  cs%u_file = -1 ; cs%v_file = -1 ; cs%cols_written = 0

end subroutine pointaccel_init

end module mom_pointaccel