MOM_PointAccel.F90

Go to the documentation of this file.

module mom_pointaccel
!***********************************************************************
!*                   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                      *
!***********************************************************************

!***********************************************************************
!*                                                                     *
!*     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.                         *
!*                                                                     *
!*  Macros written all in capital letters are defined in MOM_memory.h  *
!*                                                                     *
!*     A small fragment of the grid is shown below:                    *
!*                                                                     *
!*    j+1  x ^ x ^ x   At x:  q, CoriolisBu                            *
!*    j+1  > o > o >   At ^:  v, PFv, CAv, vh, diffv, vbt, vhtr        *
!*    j    x ^ x ^ x   At >:  u, PFu, CAu, uh, diffu, ubt, uhtr        *
!*    j    > o > o >   At o:  h, bathyT, tr, T, S                      *
!*    j-1  x ^ x ^ x                                                   *
!*        i-1  i  i+1  At x & ^:                                       *
!*           i  i+1    At > & o:                                       *
!*                                                                     *
!*  The boundaries always run through q grid points (x).               *
!*                                                                     *
!********+*********+*********+*********+*********+*********+*********+**

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_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

type, public :: pointaccel_cs ; private
  character(len=200) :: u_trunc_file ! The complete path to files in which a
  character(len=200) :: v_trunc_file ! column's worth of accelerations are
                                     ! written if velocity truncations occur.
  integer :: u_file, v_file ! The unit numbers for opened u- or v- truncation
                            ! files, or -1 if they have 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 ! A pointer to the ocean model's clock.
  type(diag_ctrl), pointer :: diag ! A pointer to a structure of shareable
                            ! ocean diagnostic fields.
! 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(), v_av => null(), & ! Time average velocities in m s-1.
    u_prev => null(), v_prev => null(), & ! Previous velocities in m s-1.
    t => null(), s => null(), &     ! Temperature and salinity in C and psu.
    pbce => null(), &               ! pbce times eta gives the baroclinic
                                    ! pressure anomaly in each layer due to
                                    ! free surface height anomalies.
                                    ! pbce has units of m s-2.
    u_accel_bt => null(), &         ! Barotropic acclerations in m s-2.
    v_accel_bt => null()

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, G, GV, CS, &
                         maxvel, minvel, 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.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                               intent(in) :: um  !< The new zonal velocity, in m s-1.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  &
                               intent(in) :: hin !< The layer thickness, in m.
  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  !< The ocean dynamics time step, in s.
  type(pointaccel_cs),         pointer    :: CS  !< The control structure returned by a previous
                                                 !! call to PointAccel_init.
  real,                        intent(in) :: maxvel, minvel
  real, optional,              intent(in) :: str !< The surface wind stress integrated over a time
                                                 !! step, in m2 s-1.
  real, dimension(SZIB_(G),SZK_(G)),         &
                     optional, intent(in) :: a   !< The layer coupling coefficients from
                                                 !! vertvisc, m.
  real, dimension(SZIB_(G),SZK_(G)),         &
                     optional, intent(in) :: hv  !< The layer thicknesses at velocity grid points,
                                                 !! from vertvisc, in m.

! 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.

! Arguments: I - The zonal index of the column to be documented.
!  (in)      j - The meridional index of the column to be documented.
!  (in)      um - The new zonal velocity, in m s-1.
!  (in)      hin - The layer thickness, in m.
!  (in)      ADp - A structure pointing to the various accelerations in
!                  the momentum equations.
!  (in)      CDp - A structure with pointers to various terms in the continuity
!                  equations.
!  (in)      dt - The model's dynamics time step.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      CS - The control structure returned by a previous call to
!                 PointAccel_init.
!  (in)      str - The surface wind stress integrated over a time
!                  step, in m2 s-1.
!  (in)      a - The layer coupling coefficients from vertvisc, m.
!  (in)      hv - The layer thicknesses at velocity grid points, from
!                 vertvisc, in m.

  real    :: f_eff, CFL
  real    :: Angstrom
  real    :: truncvel, du
  real    :: Inorm(szk_(g))
  real    :: e(szk_(g)+1)
  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 + gv%H_subroundoff

!  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)) >= maxvel) .or. &
           (min(cs%u_av(i,j,k),um(i,j,k)) <= minvel)) .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)) >= maxvel) .or. &
           (min(cs%u_av(i,j,k), um(i,j,k)) <= minvel)) .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," ",$)') (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," ",$)') (cs%u_prev(i,j,k)); enddo
    endif
    write(file,'(/,"u(3):  ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (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)) * 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)) * 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," ",$)') &
                                      ((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*adp%CAu(i,j,k)); enddo
    write(file,'(/,"PFu:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (dt*adp%PFu(i,j,k)); enddo
    write(file,'(/,"diffu: ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (dt*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*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*(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," ",$)') &
          (um(i,j,k)-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*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," ",$)') &
                                      (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,k); enddo
    endif
    if (present(hv)) then
      write(file,'(/,"hvel:  ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') hv(i,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*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," ",$)') (hin(i,j-1,k)); enddo
    write(file,'(/,"h+-:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (hin(i+1,j-1,k)); enddo
    write(file,'(/,"h-0:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (hin(i,j,k)); enddo
    write(file,'(/,"h+0:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (hin(i+1,j,k)); enddo
    write(file,'(/,"h-+:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (hin(i,j+1,k)); enddo
    write(file,'(/,"h++:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (hin(i+1,j+1,k)); enddo


    e(nz+1) = -g%bathyT(i,j)
    do k=nz,1,-1 ; e(k) = e(k+1) + 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) = -g%bathyT(i+1,j)
    do k=nz,1,-1 ; e(k) = e(k+1) + 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," ",$)') &
                                    (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*cs%v_av(i,j-1,k)*(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)*(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," ",$)') &
                                    (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*cs%v_av(i,j,k)*(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)*(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," ",$)') &
                                      (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*cs%v_av(i+1,j-1,k)*(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)*(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," ",$)') &
                          (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*cs%v_av(i+1,j,k)*(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*cs%v_av(i+1,j,k)*(hin(i+1,j,k) + hin(i+1,j+1,k))); enddo
    endif

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

  !  From here on, the normalized accelerations are written.
    if (prev_avail) then
      do k=ks,ke
        du = 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," ",$)') &
                        ((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*adp%CAu(i,j,k)*inorm(k)); enddo

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

      write(file,'(/,"diffu: ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') &
                                      (dt*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*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*(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*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," ",$)') &
            (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*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, G, GV, CS, &
                         maxvel, minvel, 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.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                               intent(in) :: vm  !< The new meridional velocity, in m s-1.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  &
                               intent(in) :: hin !< The layer thickness, in m.
  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  !< The ocean dynamics time step, in s.
  type(pointaccel_cs),         pointer    :: CS  !< The control structure returned by a previous
                                                 !! call to PointAccel_init.
  real,                        intent(in) :: maxvel, minvel
  real, optional,              intent(in) :: str !< The surface wind stress integrated over a time
                                                 !! step, in m2 s-1.
  real, dimension(SZI_(G),SZK_(G)),          &
                     optional, intent(in) :: a   !< The layer coupling coefficients from
                                                 !! vertvisc, m.
  real, dimension(SZI_(G),SZK_(G)),          &
                     optional, intent(in) :: hv  !< The layer thicknesses at velocity grid points,
                                                 !! from vertvisc, in m.

! 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.

! Arguments: i - The zonal index of the column to be documented.
!  (in)      J - The meridional index of the column to be documented.
!  (in)      vm - The new meridional velocity, in m s-1.
!  (in)      hin - The layer thickness, in m.
!  (in)      ADp - A structure pointing to the various accelerations in
!                  the momentum equations.
!  (in)      CDp - A structure with pointers to various terms in the continuity
!                  equations.
!  (in)      dt - The model's dynamics time step.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      CS - The control structure returned by a previous call to
!                 PointAccel_init.
!  (in)      str - The surface wind stress integrated over a time
!                  step, in m2 s-1.
!  (in)      a - The layer coupling coefficients from vertvisc, m.
!  (in)      hv - The layer thicknesses at velocity grid points, from
!                 vertvisc, in m.

  real    :: f_eff, CFL
  real    :: Angstrom
  real    :: truncvel, dv
  real    :: Inorm(szk_(g))
  real    :: e(szk_(g)+1)
  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 + gv%H_subroundoff

!  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)) >= maxvel) .or. &
           (min(cs%v_av(i,j,k), vm(i,j,k)) <= minvel)) .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)) >= maxvel) .or. &
           (min(cs%v_av(i,j,k), vm(i,j,k)) <= minvel)) .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," ",$)') (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," ",$)') (cs%v_prev(i,j,k)); enddo
    endif

    write(file,'(/,"v(3):  ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (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)) * 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)) * 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," ",$)') &
                                      ((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*adp%CAv(i,j,k)); enddo

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

    write(file,'(/,"diffv: ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') (dt*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*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*(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," ",$)') &
          (vm(i,j,k)-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*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," ",$)') &
                                      (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,k); enddo
    endif
    if (present(hv)) then
      write(file,'(/,"hvel:  ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') hv(i,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*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," ",$)') hin(i-1,j,k); enddo
    write(file,'(/,"h0-:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') hin(i,j,k); enddo
    write(file,'(/,"h+-:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') hin(i+1,j,k); enddo
    write(file,'(/,"h-+:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') hin(i-1,j+1,k); enddo
    write(file,'(/,"h0+:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') hin(i,j+1,k); enddo
    write(file,'(/,"h++:   ",$)')
    do k=ks,ke ; if (do_k(k)) write(file,'(ES10.3," ",$)') hin(i+1,j+1,k); enddo

    e(nz+1) = -g%bathyT(i,j)
    do k=nz,1,-1 ; e(k) = e(k+1) + 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) = -g%bathyT(i,j+1)
    do k=nz,1,-1 ; e(k) = e(k+1) + 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," ",$)') &
                                    (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," ",$)') &
            (cs%u_av(i-1,j,k) * 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," ",$)') &
                          (0.5*cs%u_prev(i-1,j,k)*(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," ",$)') &
                                    (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," ",$)') &
            (cs%u_av(i-1,j+1,k) * 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," ",$)') &
                    (0.5*cs%u_prev(i-1,j+1,k)*(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," ",$)') &
                                    (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," ",$)') &
            (cs%u_av(i,j,k) * 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," ",$)') &
                            (0.5*cs%u_prev(i,j,k)*(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," ",$)') &
                                    (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," ",$)') &
            (cs%u_av(i,j+1,k) * 0.5*(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," ",$)') &
                      (0.5*cs%u_prev(i,j+1,k)*(hin(i,j+1,k) + hin(i+1,j+1,k))); enddo
    endif

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

  !  From here on, the normalized accelerations are written.
    if (prev_avail) then
      do k=ks,ke
        dv = 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," ",$)') &
                                      ((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*adp%CAv(i,j,k)*inorm(k)); enddo
      write(file,'(/,"PFv:   ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') &
                                      (dt*adp%PFv(i,j,k)*inorm(k)); enddo
      write(file,'(/,"diffv: ",$)')
      do k=ks,ke ; if (do_k(k)) write(file,'(F10.6," ",$)') &
                                      (dt*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*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*(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*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," ",$)') &
            (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*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 needs a doxygen description
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.

! Arguments: MIS - For "MOM Internal State" a set of pointers to the fields and
!                  accelerations that make up the ocean's physical state.
!  (in)      Time - The current model time.
!  (in)      G - The ocean's grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.
!  (in)      diag - A structure that is used to regulate diagnostic output.
!  (in)      dirs - A structure containing several relevant directory paths.
!  (in/out)  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, "")
  call get_param(param_file, mdl, "U_TRUNC_FILE", cs%u_trunc_file, &
                 "The absolute path to the file where the accelerations \n"//&
                 "leading to zonal velocity truncations are written. \n"//&
                 "Leave this empty for efficiency if this diagnostic is \n"//&
                 "not needed.", default="")
  call get_param(param_file, mdl, "V_TRUNC_FILE", cs%v_trunc_file, &
                 "The absolute path to the file where the accelerations \n"//&
                 "leading to meridional velocity truncations are written. \n"//&
                 "Leave this empty for efficiency if this diagnostic is \n"//&
                 "not needed.", default="")
  call get_param(param_file, mdl, "MAX_TRUNC_FILE_SIZE_PER_PE", cs%max_writes, &
                 "The maximum number of colums of truncations that any PE \n"//&
                 "will write out during a run.", default=50)

  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)
    call log_param(param_file, mdl, "output_dir/V_TRUNC_FILE", cs%v_trunc_file)
  endif
  cs%u_file = -1 ; cs%v_file = -1 ; cs%cols_written = 0

end subroutine pointaccel_init
end module mom_pointaccel