MOM_CoriolisAdv.F90

!> Accelerations due to the Coriolis force and momentum advection
module mom_coriolisadv

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

!> \author Robert Hallberg, April 1994 - June 2002

use mom_diag_mediator, only : post_data, query_averaging_enabled, diag_ctrl
use mom_diag_mediator, only : register_diag_field, safe_alloc_ptr, time_type
use mom_error_handler, only : mom_error, mom_mesg, fatal, warning
use mom_file_parser,   only : get_param, log_version, param_file_type
use mom_grid,          only : ocean_grid_type
use mom_open_boundary, only : ocean_obc_type, obc_direction_e, obc_direction_w
use mom_open_boundary, only : obc_direction_n, obc_direction_s
use mom_string_functions, only : uppercase
use mom_unit_scaling,  only : unit_scale_type
use mom_variables,     only : accel_diag_ptrs
use mom_verticalgrid,  only : verticalgrid_type

implicit none ; private

public coradcalc, coriolisadv_init, coriolisadv_end

#include <MOM_memory.h>

!> Control structure for mom_coriolisadv
type, public :: coriolisadv_cs ; private
  integer :: coriolis_scheme !< Selects the discretization for the Coriolis terms.
                             !! Valid values are:
                             !! - SADOURNY75_ENERGY - Sadourny, 1975
                             !! - ARAKAWA_HSU90     - Arakawa & Hsu, 1990, Energy & non-div. Enstrophy
                             !! - ROBUST_ENSTRO     - Pseudo-enstrophy scheme
                             !! - SADOURNY75_ENSTRO - Sadourny, JAS 1975, Enstrophy
                             !! - ARAKAWA_LAMB81    - Arakawa & Lamb, MWR 1981, Energy & Enstrophy
                             !! - ARAKAWA_LAMB_BLEND - A blend of Arakawa & Lamb with Arakawa & Hsu and Sadourny energy.
                             !! The default, SADOURNY75_ENERGY, is the safest choice then the
                             !! deformation radius is poorly resolved.
  integer :: ke_scheme       !< KE_SCHEME selects the discretization for
                             !! the kinetic energy. Valid values are:
                             !!  KE_ARAKAWA, KE_SIMPLE_GUDONOV, KE_GUDONOV
  integer :: pv_adv_scheme   !< PV_ADV_SCHEME selects the discretization for PV advection
                             !! Valid values are:
                             !! - PV_ADV_CENTERED - centered (aka Sadourny, 75)
                             !! - PV_ADV_UPWIND1  - upwind, first order
  real    :: f_eff_max_blend !< The factor by which the maximum effective Coriolis
                             !! acceleration from any point can be increased when
                             !! blending different discretizations with the
                             !! ARAKAWA_LAMB_BLEND Coriolis scheme.  This must be
                             !! greater than 2.0, and is 4.0 by default.
  real    :: wt_lin_blend    !< A weighting value beyond which the blending between
                             !! Sadourny and Arakawa & Hsu goes linearly to 0.
                             !! This must be between 1 and 1e-15, often 1/8.
  logical :: no_slip         !< If true, no slip boundary conditions are used.
                             !! Otherwise free slip boundary conditions are assumed.
                             !! The implementation of the free slip boundary
                             !! conditions on a C-grid is much cleaner than the
                             !! no slip boundary conditions. The use of free slip
                             !! b.c.s is strongly encouraged. The no slip b.c.s
                             !! are not implemented with the biharmonic viscosity.
  logical :: bound_coriolis  !< If true, the Coriolis terms at u points are
                             !! bounded by the four estimates of (f+rv)v from the
                             !! four neighboring v points, and similarly at v
                             !! points.  This option would have no effect on the
                             !! SADOURNY75_ENERGY scheme if it were possible to
                             !! use centered difference thickness fluxes.
  logical :: coriolis_en_dis !< If CORIOLIS_EN_DIS is defined, two estimates of
                             !! the thickness fluxes are used to estimate the
                             !! Coriolis term, and the one that dissipates energy
                             !! relative to the other one is used.  This is only
                             !! available at present if Coriolis scheme is
                             !! SADOURNY75_ENERGY.
  type(time_type), pointer :: time !< A pointer to the ocean model's clock.
  type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the timing of diagnostic output.
  !>@{ Diagnostic IDs
  integer :: id_rv = -1, id_pv = -1, id_gkeu = -1, id_gkev = -1
  integer :: id_rvxu = -1, id_rvxv = -1
  ! integer :: id_hf_gKEu    = -1, id_hf_gKEv    = -1
  integer :: id_hf_gkeu_2d = -1, id_hf_gkev_2d = -1
  ! integer :: id_hf_rvxu    = -1, id_hf_rvxv    = -1
  integer :: id_hf_rvxu_2d = -1, id_hf_rvxv_2d = -1
  !>@}
end type coriolisadv_cs

!>@{ Enumeration values for Coriolis_Scheme
integer, parameter :: sadourny75_energy = 1
integer, parameter :: arakawa_hsu90     = 2
integer, parameter :: robust_enstro     = 3
integer, parameter :: sadourny75_enstro = 4
integer, parameter :: arakawa_lamb81    = 5
integer, parameter :: al_blend          = 6
character*(20), parameter :: sadourny75_energy_string = "SADOURNY75_ENERGY"
character*(20), parameter :: arakawa_hsu_string = "ARAKAWA_HSU90"
character*(20), parameter :: robust_enstro_string = "ROBUST_ENSTRO"
character*(20), parameter :: sadourny75_enstro_string = "SADOURNY75_ENSTRO"
character*(20), parameter :: arakawa_lamb_string = "ARAKAWA_LAMB81"
character*(20), parameter :: al_blend_string = "ARAKAWA_LAMB_BLEND"
!>@}
!>@{ Enumeration values for KE_Scheme
integer, parameter :: ke_arakawa        = 10
integer, parameter :: ke_simple_gudonov = 11
integer, parameter :: ke_gudonov        = 12
character*(20), parameter :: ke_arakawa_string = "KE_ARAKAWA"
character*(20), parameter :: ke_simple_gudonov_string = "KE_SIMPLE_GUDONOV"
character*(20), parameter :: ke_gudonov_string = "KE_GUDONOV"
!>@}
!>@{ Enumeration values for PV_Adv_Scheme
integer, parameter :: pv_adv_centered   = 21
integer, parameter :: pv_adv_upwind1    = 22
character*(20), parameter :: pv_adv_centered_string = "PV_ADV_CENTERED"
character*(20), parameter :: pv_adv_upwind1_string = "PV_ADV_UPWIND1"
!>@}

contains

!> Calculates the Coriolis and momentum advection contributions to the acceleration.
subroutine coradcalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS)
  type(ocean_grid_type),                     intent(in)    :: g  !< Ocen grid structure
  type(verticalgrid_type),                   intent(in)    :: gv !< Vertical grid structure
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in)    :: u  !< Zonal velocity [L T-1 ~> m s-1]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in)    :: v  !< Meridional velocity [L T-1 ~> m s-1]
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h  !< Layer thickness [H ~> m or kg m-2]
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in)    :: uh !< Zonal transport u*h*dy
                                                                 !! [H L2 T-1 ~> m3 s-1 or kg s-1]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in)    :: vh !< Meridional transport v*h*dx
                                                                 !! [H L2 T-1 ~> m3 s-1 or kg s-1]
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out)   :: cau !< Zonal acceleration due to Coriolis
                                                                  !! and momentum advection [L T-2 ~> m s-2].
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out)   :: cav !< Meridional acceleration due to Coriolis
                                                                  !! and momentum advection [L T-2 ~> m s-2].
  type(ocean_obc_type),                      pointer       :: obc !< Open boundary control structure
  type(accel_diag_ptrs),                     intent(inout) :: ad  !< Storage for acceleration diagnostics
  type(unit_scale_type),                     intent(in)    :: us  !< A dimensional unit scaling type
  type(coriolisadv_cs),                      pointer       :: cs  !< Control structure for MOM_CoriolisAdv

  ! Local variables
  real, dimension(SZIB_(G),SZJB_(G)) :: &
    q, &        ! Layer potential vorticity [H-1 T-1 ~> m-1 s-1 or m2 kg-1 s-1].
    ih_q, &     ! The inverse of thickness interpolated to q points [H-1 ~> m-1 or m2 kg-1].
    area_q      ! The sum of the ocean areas at the 4 adjacent thickness points [L2 ~> m2].

  real, dimension(SZIB_(G),SZJ_(G)) :: &
    a, b, c, d  ! a, b, c, & d are combinations of the potential vorticities
                ! surrounding an h grid point.  At small scales, a = q/4,
                ! b = q/4, etc.  All are in [H-1 T-1 ~> m-1 s-1 or m2 kg-1 s-1],
                ! and use the indexing of the corresponding u point.

  real, dimension(SZI_(G),SZJ_(G)) :: &
    area_h, &   ! The ocean area at h points [L2 ~> m2].  Area_h is used to find the
                ! average thickness in the denominator of q.  0 for land points.
    ke          ! Kinetic energy per unit mass [L2 T-2 ~> m2 s-2], KE = (u^2 + v^2)/2.
  real, dimension(SZIB_(G),SZJ_(G)) :: &
    harea_u, &  ! The cell area weighted thickness interpolated to u points
                ! times the effective areas [H L2 ~> m3 or kg].
    kex, &      ! The zonal gradient of Kinetic energy per unit mass [L T-2 ~> m s-2],
                ! KEx = d/dx KE.
    uh_center   ! Transport based on arithmetic mean h at u-points [H L2 T-1 ~> m3 s-1 or kg s-1]
  real, dimension(SZI_(G),SZJB_(G)) :: &
    harea_v, &  ! The cell area weighted thickness interpolated to v points
                ! times the effective areas [H L2 ~> m3 or kg].
    key, &      ! The meridonal gradient of Kinetic energy per unit mass [L T-2 ~> m s-2],
                ! KEy = d/dy KE.
    vh_center   ! Transport based on arithmetic mean h at v-points [H L2 T-1 ~> m3 s-1 or kg s-1]
  real, dimension(SZI_(G),SZJ_(G)) :: &
    uh_min, uh_max, &   ! The smallest and largest estimates of the volume
    vh_min, vh_max, &   ! fluxes through the faces (i.e. u*h*dy & v*h*dx)
                        ! [H L2 T-1 ~> m3 s-1 or kg s-1].
    ep_u, ep_v  ! Additional pseudo-Coriolis terms in the Arakawa and Lamb
                ! discretization [H-1 s-1 ~> m-1 s-1 or m2 kg-1 s-1].
  real, dimension(SZIB_(G),SZJB_(G)) :: &
    dvdx, dudy, & ! Contributions to the circulation around q-points [L2 T-1 ~> m2 s-1]
    abs_vort, & ! Absolute vorticity at q-points [T-1 ~> s-1].
    q2, &       ! Relative vorticity over thickness [H-1 T-1 ~> m-1 s-1 or m2 kg-1 s-1].
    max_fvq, &  ! The maximum of the adjacent values of (-u) times absolute vorticity [L T-2 ~> m s-2].
    min_fvq, &  ! The minimum of the adjacent values of (-u) times absolute vorticity [L T-2 ~> m s-2].
    max_fuq, &  ! The maximum of the adjacent values of u times absolute vorticity [L T-2 ~> m s-2].
    min_fuq     ! The minimum of the adjacent values of u times absolute vorticity [L T-2 ~> m s-2].
  real, dimension(SZIB_(G),SZJB_(G),SZK_(G)) :: &
    pv, &       ! A diagnostic array of the potential vorticities [H-1 T-1 ~> m-1 s-1 or m2 kg-1 s-1].
    rv          ! A diagnostic array of the relative vorticities [T-1 ~> s-1].
  real :: fv1, fv2, fu1, fu2   ! (f+rv)*v or (f+rv)*u [L T-2 ~> m s-2].
  real :: max_fv, max_fu       ! The maximum or minimum of the neighboring Coriolis
  real :: min_fv, min_fu       ! accelerations [L T-2 ~> m s-2], i.e. max(min)_fu(v)q.

  real, parameter :: c1_12=1.0/12.0 ! C1_12 = 1/12
  real, parameter :: c1_24=1.0/24.0 ! C1_24 = 1/24
  real :: absolute_vorticity     ! Absolute vorticity [T-1 ~> s-1].
  real :: relative_vorticity     ! Relative vorticity [T-1 ~> s-1].
  real :: ih                     ! Inverse of thickness [H-1 ~> m-1 or m2 kg-1].
  real :: max_ihq, min_ihq       ! The maximum and minimum of the nearby Ihq [H-1 ~> m-1 or m2 kg-1].
  real :: harea_q                ! The sum of area times thickness of the cells
                                 ! surrounding a q point [H L2 ~> m3 or kg].
  real :: h_neglect              ! A thickness that is so small it is usually
                                 ! lost in roundoff and can be neglected [H ~> m or kg m-2].
  real :: temp1, temp2           ! Temporary variables [L2 T-2 ~> m2 s-2].
  real :: eps_vel                ! A tiny, positive velocity [L T-1 ~> m s-1].

  real :: uhc, vhc               ! Centered estimates of uh and vh [H L2 T-1 ~> m3 s-1 or kg s-1].
  real :: uhm, vhm               ! The input estimates of uh and vh [H L2 T-1 ~> m3 s-1 or kg s-1].
  real :: c1, c2, c3, slope      ! Nondimensional parameters for the Coriolis limiter scheme.

  real :: fe_m2         ! Nondimensional temporary variables asssociated with
  real :: rat_lin       ! the ARAKAWA_LAMB_BLEND scheme.
  real :: rat_m1        ! The ratio of the maximum neighboring inverse thickness
                        ! to the minimum inverse thickness minus 1. rat_m1 >= 0.
  real :: al_wt         ! The relative weight of the Arakawa & Lamb scheme to the
                        ! Arakawa & Hsu scheme, nondimensional between 0 and 1.
  real :: sad_wt        ! The relative weight of the Sadourny energy scheme to
                        ! the other two with the ARAKAWA_LAMB_BLEND scheme,
                        ! nondimensional between 0 and 1.

  real :: heff1, heff2  ! Temporary effective H at U or V points [H ~> m or kg m-2].
  real :: heff3, heff4  ! Temporary effective H at U or V points [H ~> m or kg m-2].
  real :: h_tiny        ! A very small thickness [H ~> m or kg m-2].
  real :: uheff, vheff  ! More temporary variables [H L2 T-1 ~> m3 s-1 or kg s-1].
  real :: quheff,qvheff ! More temporary variables [H L2 T-1 s-1 ~> m3 s-2 or kg s-2].
  integer :: i, j, k, n, is, ie, js, je, isq, ieq, jsq, jeq, nz

! Diagnostics for fractional thickness-weighted terms
  real, allocatable, dimension(:,:) :: &
    hf_gkeu_2d, hf_gkev_2d, & ! Depth sum of hf_gKEu, hf_gKEv [L T-2 ~> m s-2].
    hf_rvxu_2d, hf_rvxv_2d    ! Depth sum of hf_rvxu, hf_rvxv [L T-2 ~> m s-2].
  !real, allocatable, dimension(:,:,:) :: &
  !  hf_gKEu, hf_gKEv, & ! accel. due to KE gradient x fract. thickness  [L T-2 ~> m s-2].
  !  hf_rvxu, hf_rvxv    ! accel. due to RV x fract. thickness [L T-2 ~> m s-2].
  ! 3D diagnostics hf_gKEu etc. are commented because there is no clarity on proper remapping grid option.
  ! The code is retained for degugging purposes in the future.

! To work, the following fields must be set outside of the usual
! is to ie range before this subroutine is called:
!   v(is-1:ie+2,js-1:je+1), u(is-1:ie+1,js-1:je+2), h(is-1:ie+2,js-1:je+2),
!   uh(is-1,ie,js:je+1) and vh(is:ie+1,js-1:je).

  if (.not.associated(cs)) call mom_error(fatal, &
         "MOM_CoriolisAdv: Module must be initialized before it is used.")
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB ; nz = g%ke
  h_neglect = gv%H_subroundoff
  eps_vel = 1.0e-10*us%m_s_to_L_T
  h_tiny = gv%Angstrom_H  ! Perhaps this should be set to h_neglect instead.

  !$OMP parallel do default(private) shared(Isq,Ieq,Jsq,Jeq,G,Area_h)
  do j=jsq-1,jeq+2 ; do i=isq-1,ieq+2
    area_h(i,j) = g%mask2dT(i,j) * g%areaT(i,j)
  enddo ; enddo
  if (associated(obc)) then ; do n=1,obc%number_of_segments
    if (.not. obc%segment(n)%on_pe) cycle
    i = obc%segment(n)%HI%IsdB ; j = obc%segment(n)%HI%JsdB
    if (obc%segment(n)%is_N_or_S .and. (j >= jsq-1) .and. (j <= jeq+1)) then
      do i = max(isq-1,obc%segment(n)%HI%isd), min(ieq+2,obc%segment(n)%HI%ied)
        if (obc%segment(n)%direction == obc_direction_n) then
          area_h(i,j+1) = area_h(i,j)
        else ! (OBC%segment(n)%direction == OBC_DIRECTION_S)
          area_h(i,j) = area_h(i,j+1)
        endif
      enddo
    elseif (obc%segment(n)%is_E_or_W .and. (i >= isq-1) .and. (i <= ieq+1)) then
      do j = max(jsq-1,obc%segment(n)%HI%jsd), min(jeq+2,obc%segment(n)%HI%jed)
        if (obc%segment(n)%direction == obc_direction_e) then
          area_h(i+1,j) = area_h(i,j)
        else ! (OBC%segment(n)%direction == OBC_DIRECTION_W)
          area_h(i,j) = area_h(i+1,j)
        endif
      enddo
    endif
  enddo ; endif
  !$OMP parallel do default(private) shared(Isq,Ieq,Jsq,Jeq,G,Area_h,Area_q)
  do j=jsq-1,jeq+1 ; do i=isq-1,ieq+1
    area_q(i,j) = (area_h(i,j) + area_h(i+1,j+1)) + &
                  (area_h(i+1,j) + area_h(i,j+1))
  enddo ; enddo

  !$OMP parallel do default(private) shared(u,v,h,uh,vh,CAu,CAv,G,CS,AD,Area_h,Area_q,&
  !$OMP                        RV,PV,is,ie,js,je,Isq,Ieq,Jsq,Jeq,nz,h_neglect,h_tiny,OBC,eps_vel)
  do k=1,nz

    ! Here the second order accurate layer potential vorticities, q,
    ! are calculated.  hq is  second order accurate in space.  Relative
    ! vorticity is second order accurate everywhere with free slip b.c.s,
    ! but only first order accurate at boundaries with no slip b.c.s.
    ! First calculate the contributions to the circulation around the q-point.
    do j=jsq-1,jeq+1 ; do i=isq-1,ieq+1
      dvdx(i,j) = (v(i+1,j,k)*g%dyCv(i+1,j) - v(i,j,k)*g%dyCv(i,j))
      dudy(i,j) = (u(i,j+1,k)*g%dxCu(i,j+1) - u(i,j,k)*g%dxCu(i,j))
    enddo ; enddo
    do j=jsq-1,jeq+1 ; do i=isq-1,ieq+2
      harea_v(i,j) = 0.5*(area_h(i,j) * h(i,j,k) + area_h(i,j+1) * h(i,j+1,k))
    enddo ; enddo
    do j=jsq-1,jeq+2 ; do i=isq-1,ieq+1
      harea_u(i,j) = 0.5*(area_h(i,j) * h(i,j,k) + area_h(i+1,j) * h(i+1,j,k))
    enddo ; enddo
    if (cs%Coriolis_En_Dis) then
      do j=jsq,jeq+1 ; do i=is-1,ie
        uh_center(i,j) = 0.5 * (g%dy_Cu(i,j) * u(i,j,k)) * (h(i,j,k) + h(i+1,j,k))
      enddo ; enddo
      do j=js-1,je ; do i=isq,ieq+1
        vh_center(i,j) = 0.5 * (g%dx_Cv(i,j) * v(i,j,k)) * (h(i,j,k) + h(i,j+1,k))
      enddo ; enddo
    endif

    ! Adjust circulation components to relative vorticity and thickness projected onto
    ! velocity points on open boundaries.
    if (associated(obc)) then ; do n=1,obc%number_of_segments
      if (.not. obc%segment(n)%on_pe) cycle
      i = obc%segment(n)%HI%IsdB ; j = obc%segment(n)%HI%JsdB
      if (obc%segment(n)%is_N_or_S .and. (j >= jsq-1) .and. (j <= jeq+1)) then
        if (obc%zero_vorticity) then ; do i=obc%segment(n)%HI%IsdB,obc%segment(n)%HI%IedB
          dvdx(i,j) = 0. ; dudy(i,j) = 0.
        enddo ; endif
        if (obc%freeslip_vorticity) then ; do i=obc%segment(n)%HI%IsdB,obc%segment(n)%HI%IedB
          dudy(i,j) = 0.
        enddo ; endif
        if (obc%computed_vorticity) then ; do i=obc%segment(n)%HI%IsdB,obc%segment(n)%HI%IedB
          if (obc%segment(n)%direction == obc_direction_n) then
            dudy(i,j) = 2.0*(obc%segment(n)%tangential_vel(i,j,k) - u(i,j,k))*g%dxCu(i,j)
          else ! (OBC%segment(n)%direction == OBC_DIRECTION_S)
            dudy(i,j) = 2.0*(u(i,j+1,k) - obc%segment(n)%tangential_vel(i,j,k))*g%dxCu(i,j+1)
          endif
        enddo ; endif
        if (obc%specified_vorticity) then ; do i=obc%segment(n)%HI%IsdB,obc%segment(n)%HI%IedB
          if (obc%segment(n)%direction == obc_direction_n) then
            dudy(i,j) = obc%segment(n)%tangential_grad(i,j,k)*g%dxCu(i,j)*g%dyBu(i,j)
          else ! (OBC%segment(n)%direction == OBC_DIRECTION_S)
            dudy(i,j) = obc%segment(n)%tangential_grad(i,j,k)*g%dxCu(i,j+1)*g%dyBu(i,j)
          endif
        enddo ; endif

        ! Project thicknesses across OBC points with a no-gradient condition.
        do i = max(isq-1,obc%segment(n)%HI%isd), min(ieq+2,obc%segment(n)%HI%ied)
          if (obc%segment(n)%direction == obc_direction_n) then
            harea_v(i,j) = 0.5 * (area_h(i,j) + area_h(i,j+1)) * h(i,j,k)
          else ! (OBC%segment(n)%direction == OBC_DIRECTION_S)
            harea_v(i,j) = 0.5 * (area_h(i,j) + area_h(i,j+1)) * h(i,j+1,k)
          endif
        enddo

        if (cs%Coriolis_En_Dis) then
          do i = max(isq-1,obc%segment(n)%HI%isd), min(ieq+2,obc%segment(n)%HI%ied)
            if (obc%segment(n)%direction == obc_direction_n) then
              vh_center(i,j) = g%dx_Cv(i,j) * v(i,j,k) * h(i,j,k)
            else ! (OBC%segment(n)%direction == OBC_DIRECTION_S)
              vh_center(i,j) = g%dx_Cv(i,j) * v(i,j,k) * h(i,j+1,k)
            endif
          enddo
        endif
      elseif (obc%segment(n)%is_E_or_W .and. (i >= isq-1) .and. (i <= ieq+1)) then
        if (obc%zero_vorticity) then ; do j=obc%segment(n)%HI%JsdB,obc%segment(n)%HI%JedB
          dvdx(i,j) = 0. ; dudy(i,j) = 0.
        enddo ; endif
        if (obc%freeslip_vorticity) then ; do j=obc%segment(n)%HI%JsdB,obc%segment(n)%HI%JedB
          dvdx(i,j) = 0.
        enddo ; endif
        if (obc%computed_vorticity) then ; do j=obc%segment(n)%HI%JsdB,obc%segment(n)%HI%JedB
          if (obc%segment(n)%direction == obc_direction_e) then
            dvdx(i,j) = 2.0*(obc%segment(n)%tangential_vel(i,j,k) - v(i,j,k))*g%dyCv(i,j)
          else ! (OBC%segment(n)%direction == OBC_DIRECTION_W)
            dvdx(i,j) = 2.0*(v(i+1,j,k) - obc%segment(n)%tangential_vel(i,j,k))*g%dyCv(i+1,j)
          endif
        enddo ; endif
        if (obc%specified_vorticity) then ; do j=obc%segment(n)%HI%JsdB,obc%segment(n)%HI%JedB
          if (obc%segment(n)%direction == obc_direction_e) then
            dvdx(i,j) = obc%segment(n)%tangential_grad(i,j,k)*g%dyCv(i,j)*g%dxBu(i,j)
          else ! (OBC%segment(n)%direction == OBC_DIRECTION_W)
            dvdx(i,j) = obc%segment(n)%tangential_grad(i,j,k)*g%dyCv(i+1,j)*g%dxBu(i,j)
          endif
        enddo ; endif

        ! Project thicknesses across OBC points with a no-gradient condition.
        do j = max(jsq-1,obc%segment(n)%HI%jsd), min(jeq+2,obc%segment(n)%HI%jed)
          if (obc%segment(n)%direction == obc_direction_e) then
            harea_u(i,j) = 0.5*(area_h(i,j) + area_h(i+1,j)) * h(i,j,k)
          else ! (OBC%segment(n)%direction == OBC_DIRECTION_W)
            harea_u(i,j) = 0.5*(area_h(i,j) + area_h(i+1,j)) * h(i+1,j,k)
          endif
        enddo
        if (cs%Coriolis_En_Dis) then
          do j = max(jsq-1,obc%segment(n)%HI%jsd), min(jeq+2,obc%segment(n)%HI%jed)
            if (obc%segment(n)%direction == obc_direction_e) then
              uh_center(i,j) = g%dy_Cu(i,j) * u(i,j,k) * h(i,j,k)
            else ! (OBC%segment(n)%direction == OBC_DIRECTION_W)
              uh_center(i,j) = g%dy_Cu(i,j) * u(i,j,k) * h(i+1,j,k)
            endif
          enddo
        endif
      endif
    enddo ; endif

    if (associated(obc)) then ; do n=1,obc%number_of_segments
      if (.not. obc%segment(n)%on_pe) cycle
      ! Now project thicknesses across cell-corner points in the OBCs.  The two
      ! projections have to occur in sequence and can not be combined easily.
      i = obc%segment(n)%HI%IsdB ; j = obc%segment(n)%HI%JsdB
      if (obc%segment(n)%is_N_or_S .and. (j >= jsq-1) .and. (j <= jeq+1)) then
        do i = max(isq-1,obc%segment(n)%HI%IsdB), min(ieq+1,obc%segment(n)%HI%IedB)
          if (obc%segment(n)%direction == obc_direction_n) then
            if (area_h(i,j) + area_h(i+1,j) > 0.0) then
              harea_u(i,j+1) = harea_u(i,j) * ((area_h(i,j+1) + area_h(i+1,j+1)) / &
                                               (area_h(i,j) + area_h(i+1,j)))
            else ; harea_u(i,j+1) = 0.0 ; endif
          else ! (OBC%segment(n)%direction == OBC_DIRECTION_S)
            if (area_h(i,j+1) + area_h(i+1,j+1) > 0.0) then
              harea_u(i,j) = harea_u(i,j+1) * ((area_h(i,j) + area_h(i+1,j)) / &
                                               (area_h(i,j+1) + area_h(i+1,j+1)))
            else ; harea_u(i,j) = 0.0 ; endif
          endif
        enddo
      elseif (obc%segment(n)%is_E_or_W .and. (i >= isq-1) .and. (i <= ieq+1)) then
        do j = max(jsq-1,obc%segment(n)%HI%JsdB), min(jeq+1,obc%segment(n)%HI%JedB)
          if (obc%segment(n)%direction == obc_direction_e) then
            if (area_h(i,j) + area_h(i,j+1) > 0.0) then
              harea_v(i+1,j) = harea_v(i,j) * ((area_h(i+1,j) + area_h(i+1,j+1)) / &
                                               (area_h(i,j) + area_h(i,j+1)))
            else ; harea_v(i+1,j) = 0.0 ; endif
          else ! (OBC%segment(n)%direction == OBC_DIRECTION_W)
            harea_v(i,j) = 0.5 * (area_h(i,j) + area_h(i,j+1)) * h(i,j+1,k)
            if (area_h(i+1,j) + area_h(i+1,j+1) > 0.0) then
              harea_v(i,j) = harea_v(i+1,j) * ((area_h(i,j) + area_h(i,j+1)) / &
                                               (area_h(i+1,j) + area_h(i+1,j+1)))
            else ; harea_v(i,j) = 0.0 ; endif
          endif
        enddo
      endif
    enddo ; endif

    do j=jsq-1,jeq+1 ; do i=isq-1,ieq+1
      if (cs%no_slip ) then
        relative_vorticity = (2.0-g%mask2dBu(i,j)) * (dvdx(i,j) - dudy(i,j)) * g%IareaBu(i,j)
      else
        relative_vorticity = g%mask2dBu(i,j) * (dvdx(i,j) - dudy(i,j)) * g%IareaBu(i,j)
      endif
      absolute_vorticity = g%CoriolisBu(i,j) + relative_vorticity
      ih = 0.0
      if (area_q(i,j) > 0.0) then
        harea_q = (harea_u(i,j) + harea_u(i,j+1)) + (harea_v(i,j) + harea_v(i+1,j))
        ih = area_q(i,j) / (harea_q + h_neglect*area_q(i,j))
      endif
      q(i,j) = absolute_vorticity * ih
      abs_vort(i,j) = absolute_vorticity
      ih_q(i,j) = ih

      if (cs%bound_Coriolis) then
        fv1 = absolute_vorticity * v(i+1,j,k)
        fv2 = absolute_vorticity * v(i,j,k)
        fu1 = -absolute_vorticity * u(i,j+1,k)
        fu2 = -absolute_vorticity * u(i,j,k)
        if (fv1 > fv2) then
          max_fvq(i,j) = fv1 ; min_fvq(i,j) = fv2
        else
          max_fvq(i,j) = fv2 ; min_fvq(i,j) = fv1
        endif
        if (fu1 > fu2) then
          max_fuq(i,j) = fu1 ; min_fuq(i,j) = fu2
        else
          max_fuq(i,j) = fu2 ; min_fuq(i,j) = fu1
        endif
      endif

      if (cs%id_rv > 0) rv(i,j,k) = relative_vorticity
      if (cs%id_PV > 0) pv(i,j,k) = q(i,j)
      if (associated(ad%rv_x_v) .or. associated(ad%rv_x_u)) &
        q2(i,j) = relative_vorticity * ih
    enddo ; enddo

    !   a, b, c, and d are combinations of neighboring potential
    ! vorticities which form the Arakawa and Hsu vorticity advection
    ! scheme.  All are defined at u grid points.

    if (cs%Coriolis_Scheme == arakawa_hsu90) then
      do j=jsq,jeq+1
        do i=is-1,ieq
          a(i,j) = (q(i,j) + (q(i+1,j) + q(i,j-1))) * c1_12
          d(i,j) = ((q(i,j) + q(i+1,j-1)) + q(i,j-1)) * c1_12
        enddo
        do i=isq,ieq
          b(i,j) = (q(i,j) + (q(i-1,j) + q(i,j-1))) * c1_12
          c(i,j) = ((q(i,j) + q(i-1,j-1)) + q(i,j-1)) * c1_12
        enddo
      enddo
    elseif (cs%Coriolis_Scheme == arakawa_lamb81) then
      do j=jsq,jeq+1 ; do i=isq,ieq+1
        a(i-1,j) = (2.0*(q(i,j) + q(i-1,j-1)) + (q(i-1,j) + q(i,j-1))) * c1_24
        d(i-1,j) = ((q(i,j) + q(i-1,j-1)) + 2.0*(q(i-1,j) + q(i,j-1))) * c1_24
        b(i,j) =   ((q(i,j) + q(i-1,j-1)) + 2.0*(q(i-1,j) + q(i,j-1))) * c1_24
        c(i,j) =   (2.0*(q(i,j) + q(i-1,j-1)) + (q(i-1,j) + q(i,j-1))) * c1_24
        ep_u(i,j) = ((q(i,j) - q(i-1,j-1)) + (q(i-1,j) - q(i,j-1))) * c1_24
        ep_v(i,j) = (-(q(i,j) - q(i-1,j-1)) + (q(i-1,j) - q(i,j-1))) * c1_24
      enddo ; enddo
    elseif (cs%Coriolis_Scheme == al_blend) then
      fe_m2 = cs%F_eff_max_blend - 2.0
      rat_lin = 1.5 * fe_m2 / max(cs%wt_lin_blend, 1.0e-16)

      ! This allows the code to always give Sadourny Energy
      if (cs%F_eff_max_blend <= 2.0) then ; fe_m2 = -1. ; rat_lin = -1.0 ; endif

      do j=jsq,jeq+1 ; do i=isq,ieq+1
        min_ihq = min(ih_q(i-1,j-1), ih_q(i,j-1), ih_q(i-1,j), ih_q(i,j))
        max_ihq = max(ih_q(i-1,j-1), ih_q(i,j-1), ih_q(i-1,j), ih_q(i,j))
        rat_m1 = 1.0e15
        if (max_ihq < 1.0e15*min_ihq) rat_m1 = max_ihq / min_ihq - 1.0
        ! The weights used here are designed to keep the effective Coriolis
        ! acceleration from any one point on its neighbors within a factor
        ! of F_eff_max.  The minimum permitted value is 2 (the factor for
        ! Sadourny's energy conserving scheme).

        ! Determine the relative weights of Arakawa & Lamb vs. Arakawa and Hsu.
        if (rat_m1 <= fe_m2) then ; al_wt = 1.0
        elseif (rat_m1 < 1.5*fe_m2) then ; al_wt = 3.0*fe_m2 / rat_m1 - 2.0
        else ; al_wt = 0.0 ; endif

        ! Determine the relative weights of Sadourny Energy vs. the other two.
        if (rat_m1 <= 1.5*fe_m2) then ; sad_wt = 0.0
        elseif (rat_m1 <= rat_lin) then
          sad_wt = 1.0 - (1.5*fe_m2) / rat_m1
        elseif (rat_m1 < 2.0*rat_lin) then
          sad_wt = 1.0 - (cs%wt_lin_blend / rat_lin) * (rat_m1 - 2.0*rat_lin)
        else ; sad_wt = 1.0 ; endif

        a(i-1,j) = sad_wt * 0.25 * q(i-1,j) + (1.0 - sad_wt) * &
                   ( ((2.0-al_wt)* q(i-1,j) + al_wt*q(i,j-1)) + &
                      2.0 * (q(i,j) + q(i-1,j-1)) ) * c1_24
        d(i-1,j) = sad_wt * 0.25 * q(i-1,j-1) + (1.0 - sad_wt) * &
                   ( ((2.0-al_wt)* q(i-1,j-1) + al_wt*q(i,j)) + &
                      2.0 * (q(i-1,j) + q(i,j-1)) ) * c1_24
        b(i,j) =   sad_wt * 0.25 * q(i,j) + (1.0 - sad_wt) * &
                   ( ((2.0-al_wt)* q(i,j) + al_wt*q(i-1,j-1)) + &
                      2.0 * (q(i-1,j) + q(i,j-1)) ) * c1_24
        c(i,j) =   sad_wt * 0.25 * q(i,j-1) + (1.0 - sad_wt) * &
                   ( ((2.0-al_wt)* q(i,j-1) + al_wt*q(i-1,j)) + &
                      2.0 * (q(i,j) + q(i-1,j-1)) ) * c1_24
        ep_u(i,j) = al_wt  * ((q(i,j) - q(i-1,j-1)) + (q(i-1,j) - q(i,j-1))) * c1_24
        ep_v(i,j) = al_wt * (-(q(i,j) - q(i-1,j-1)) + (q(i-1,j) - q(i,j-1))) * c1_24
      enddo ; enddo
    endif

    if (cs%Coriolis_En_Dis) then
    !  c1 = 1.0-1.5*RANGE ; c2 = 1.0-RANGE ; c3 = 2.0 ; slope = 0.5
      c1 = 1.0-1.5*0.5 ; c2 = 1.0-0.5 ; c3 = 2.0 ; slope = 0.5

      do j=jsq,jeq+1 ; do i=is-1,ie
        uhc = uh_center(i,j)
        uhm = uh(i,j,k)
        ! This sometimes matters with some types of open boundary conditions.
        if (g%dy_Cu(i,j) == 0.0) uhc = uhm

        if (abs(uhc) < 0.1*abs(uhm)) then
          uhm = 10.0*uhc
        elseif (abs(uhc) > c1*abs(uhm)) then
          if (abs(uhc) < c2*abs(uhm)) then ; uhc = (3.0*uhc+(1.0-c2*3.0)*uhm)
          elseif (abs(uhc) <= c3*abs(uhm)) then ; uhc = uhm
          else ; uhc = slope*uhc+(1.0-c3*slope)*uhm
          endif
        endif

        if (uhc > uhm) then
          uh_min(i,j) = uhm ; uh_max(i,j) = uhc
        else
          uh_max(i,j) = uhm ; uh_min(i,j) = uhc
        endif
      enddo ; enddo
      do j=js-1,je ; do i=isq,ieq+1
        vhc = vh_center(i,j)
        vhm = vh(i,j,k)
        ! This sometimes matters with some types of open boundary conditions.
        if (g%dx_Cv(i,j) == 0.0) vhc = vhm

        if (abs(vhc) < 0.1*abs(vhm)) then
          vhm = 10.0*vhc
        elseif (abs(vhc) > c1*abs(vhm)) then
          if (abs(vhc) < c2*abs(vhm)) then ; vhc = (3.0*vhc+(1.0-c2*3.0)*vhm)
          elseif (abs(vhc) <= c3*abs(vhm)) then ; vhc = vhm
          else ; vhc = slope*vhc+(1.0-c3*slope)*vhm
          endif
        endif

        if (vhc > vhm) then
          vh_min(i,j) = vhm ; vh_max(i,j) = vhc
        else
          vh_max(i,j) = vhm ; vh_min(i,j) = vhc
        endif
      enddo ; enddo
    endif

    ! Calculate KE and the gradient of KE
    call gradke(u, v, h, ke, kex, key, k, obc, g, us, cs)

    ! Calculate the tendencies of zonal velocity due to the Coriolis
    ! force and momentum advection.  On a Cartesian grid, this is
    !     CAu =  q * vh - d(KE)/dx.
    if (cs%Coriolis_Scheme == sadourny75_energy) then
      if (cs%Coriolis_En_Dis) then
        ! Energy dissipating biased scheme, Hallberg 200x
        do j=js,je ; do i=isq,ieq
          if (q(i,j)*u(i,j,k) == 0.0) then
            temp1 = q(i,j) * ( (vh_max(i,j)+vh_max(i+1,j)) &
                             + (vh_min(i,j)+vh_min(i+1,j)) )*0.5
          elseif (q(i,j)*u(i,j,k) < 0.0) then
            temp1 = q(i,j) * (vh_max(i,j)+vh_max(i+1,j))
          else
            temp1 = q(i,j) * (vh_min(i,j)+vh_min(i+1,j))
          endif
          if (q(i,j-1)*u(i,j,k) == 0.0) then
            temp2 = q(i,j-1) * ( (vh_max(i,j-1)+vh_max(i+1,j-1)) &
                               + (vh_min(i,j-1)+vh_min(i+1,j-1)) )*0.5
          elseif (q(i,j-1)*u(i,j,k) < 0.0) then
            temp2 = q(i,j-1) * (vh_max(i,j-1)+vh_max(i+1,j-1))
          else
            temp2 = q(i,j-1) * (vh_min(i,j-1)+vh_min(i+1,j-1))
          endif
          cau(i,j,k) = 0.25 * g%IdxCu(i,j) * (temp1 + temp2)
        enddo ; enddo
      else
        ! Energy conserving scheme, Sadourny 1975
        do j=js,je ; do i=isq,ieq
          cau(i,j,k) = 0.25 * &
            (q(i,j) * (vh(i+1,j,k) + vh(i,j,k)) + &
             q(i,j-1) * (vh(i,j-1,k) + vh(i+1,j-1,k))) * g%IdxCu(i,j)
        enddo ; enddo
      endif
    elseif (cs%Coriolis_Scheme == sadourny75_enstro) then
      do j=js,je ; do i=isq,ieq
        cau(i,j,k) = 0.125 * (g%IdxCu(i,j) * (q(i,j) + q(i,j-1))) * &
                     ((vh(i+1,j,k) + vh(i,j,k)) + (vh(i,j-1,k) + vh(i+1,j-1,k)))
      enddo ; enddo
    elseif ((cs%Coriolis_Scheme == arakawa_hsu90) .or. &
            (cs%Coriolis_Scheme == arakawa_lamb81) .or. &
            (cs%Coriolis_Scheme == al_blend)) then
      ! (Global) Energy and (Local) Enstrophy conserving, Arakawa & Hsu 1990
      do j=js,je ; do i=isq,ieq
        cau(i,j,k) = ((a(i,j) * vh(i+1,j,k) +  c(i,j) * vh(i,j-1,k))  + &
                      (b(i,j) * vh(i,j,k) +  d(i,j) * vh(i+1,j-1,k))) * g%IdxCu(i,j)
      enddo ; enddo
    elseif (cs%Coriolis_Scheme == robust_enstro) then
      ! An enstrophy conserving scheme robust to vanishing layers
      ! Note: Heffs are in lieu of h_at_v that should be returned by the
      !       continuity solver. AJA
      do j=js,je ; do i=isq,ieq
        heff1 = abs(vh(i,j,k) * g%IdxCv(i,j)) / (eps_vel+abs(v(i,j,k)))
        heff1 = max(heff1, min(h(i,j,k),h(i,j+1,k)))
        heff1 = min(heff1, max(h(i,j,k),h(i,j+1,k)))
        heff2 = abs(vh(i,j-1,k) * g%IdxCv(i,j-1)) / (eps_vel+abs(v(i,j-1,k)))
        heff2 = max(heff2, min(h(i,j-1,k),h(i,j,k)))
        heff2 = min(heff2, max(h(i,j-1,k),h(i,j,k)))
        heff3 = abs(vh(i+1,j,k) * g%IdxCv(i+1,j)) / (eps_vel+abs(v(i+1,j,k)))
        heff3 = max(heff3, min(h(i+1,j,k),h(i+1,j+1,k)))
        heff3 = min(heff3, max(h(i+1,j,k),h(i+1,j+1,k)))
        heff4 = abs(vh(i+1,j-1,k) * g%IdxCv(i+1,j-1)) / (eps_vel+abs(v(i+1,j-1,k)))
        heff4 = max(heff4, min(h(i+1,j-1,k),h(i+1,j,k)))
        heff4 = min(heff4, max(h(i+1,j-1,k),h(i+1,j,k)))
        if (cs%PV_Adv_Scheme == pv_adv_centered) then
          cau(i,j,k) = 0.5*(abs_vort(i,j)+abs_vort(i,j-1)) * &
                       ((vh(i,j,k) + vh(i+1,j-1,k)) + (vh(i,j-1,k) + vh(i+1,j,k)) ) /  &
                       (h_tiny + ((heff1+heff4) + (heff2+heff3)) ) * g%IdxCu(i,j)
        elseif (cs%PV_Adv_Scheme == pv_adv_upwind1) then
          vheff = ((vh(i,j,k) + vh(i+1,j-1,k)) + (vh(i,j-1,k) + vh(i+1,j,k)) )
          qvheff = 0.5*( (abs_vort(i,j)+abs_vort(i,j-1))*vheff &
                        -(abs_vort(i,j)-abs_vort(i,j-1))*abs(vheff) )
          cau(i,j,k) = (qvheff / ( h_tiny + ((heff1+heff4) + (heff2+heff3)) ) ) * g%IdxCu(i,j)
        endif
      enddo ; enddo
    endif
    ! Add in the additonal terms with Arakawa & Lamb.
    if ((cs%Coriolis_Scheme == arakawa_lamb81) .or. &
        (cs%Coriolis_Scheme == al_blend)) then ; do j=js,je ; do i=isq,ieq
      cau(i,j,k) = cau(i,j,k) + &
            (ep_u(i,j)*uh(i-1,j,k) - ep_u(i+1,j)*uh(i+1,j,k)) * g%IdxCu(i,j)
    enddo ; enddo ; endif


    if (cs%bound_Coriolis) then
      do j=js,je ; do i=isq,ieq
        max_fv = max(max_fvq(i,j), max_fvq(i,j-1))
        min_fv = min(min_fvq(i,j), min_fvq(i,j-1))
       ! CAu(I,j,k) = min( CAu(I,j,k), max_fv )
       ! CAu(I,j,k) = max( CAu(I,j,k), min_fv )
        if (cau(i,j,k) > max_fv) then
            cau(i,j,k) = max_fv
        else
          if (cau(i,j,k) < min_fv) cau(i,j,k) = min_fv
        endif
      enddo ; enddo
    endif

    ! Term - d(KE)/dx.
    do j=js,je ; do i=isq,ieq
      cau(i,j,k) = cau(i,j,k) - kex(i,j)
      if (associated(ad%gradKEu)) ad%gradKEu(i,j,k) = -kex(i,j)
    enddo ; enddo


    ! Calculate the tendencies of meridional velocity due to the Coriolis
    ! force and momentum advection.  On a Cartesian grid, this is
    !     CAv = - q * uh - d(KE)/dy.
    if (cs%Coriolis_Scheme == sadourny75_energy) then
      if (cs%Coriolis_En_Dis) then
        ! Energy dissipating biased scheme, Hallberg 200x
        do j=jsq,jeq ; do i=is,ie
          if (q(i-1,j)*v(i,j,k) == 0.0) then
            temp1 = q(i-1,j) * ( (uh_max(i-1,j)+uh_max(i-1,j+1)) &
                               + (uh_min(i-1,j)+uh_min(i-1,j+1)) )*0.5
          elseif (q(i-1,j)*v(i,j,k) > 0.0) then
            temp1 = q(i-1,j) * (uh_max(i-1,j)+uh_max(i-1,j+1))
          else
            temp1 = q(i-1,j) * (uh_min(i-1,j)+uh_min(i-1,j+1))
          endif
          if (q(i,j)*v(i,j,k) == 0.0) then
            temp2 = q(i,j) * ( (uh_max(i,j)+uh_max(i,j+1)) &
                             + (uh_min(i,j)+uh_min(i,j+1)) )*0.5
          elseif (q(i,j)*v(i,j,k) > 0.0) then
            temp2 = q(i,j) * (uh_max(i,j)+uh_max(i,j+1))
          else
            temp2 = q(i,j) * (uh_min(i,j)+uh_min(i,j+1))
          endif
          cav(i,j,k) = -0.25 * g%IdyCv(i,j) * (temp1 + temp2)
        enddo ; enddo
      else
        ! Energy conserving scheme, Sadourny 1975
        do j=jsq,jeq ; do i=is,ie
          cav(i,j,k) = - 0.25* &
              (q(i-1,j)*(uh(i-1,j,k) + uh(i-1,j+1,k)) + &
               q(i,j)*(uh(i,j,k) + uh(i,j+1,k))) * g%IdyCv(i,j)
        enddo ; enddo
      endif
    elseif (cs%Coriolis_Scheme == sadourny75_enstro) then
      do j=jsq,jeq ; do i=is,ie
        cav(i,j,k) = -0.125 * (g%IdyCv(i,j) * (q(i-1,j) + q(i,j))) * &
                     ((uh(i-1,j,k) + uh(i-1,j+1,k)) + (uh(i,j,k) + uh(i,j+1,k)))
      enddo ; enddo
    elseif ((cs%Coriolis_Scheme == arakawa_hsu90) .or. &
            (cs%Coriolis_Scheme == arakawa_lamb81) .or. &
            (cs%Coriolis_Scheme == al_blend)) then
      ! (Global) Energy and (Local) Enstrophy conserving, Arakawa & Hsu 1990
      do j=jsq,jeq ; do i=is,ie
        cav(i,j,k) = - ((a(i-1,j)   * uh(i-1,j,k) + &
                         c(i,j+1)   * uh(i,j+1,k))  &
                      + (b(i,j)     * uh(i,j,k) +   &
                         d(i-1,j+1) * uh(i-1,j+1,k))) * g%IdyCv(i,j)
      enddo ; enddo
    elseif (cs%Coriolis_Scheme == robust_enstro) then
      ! An enstrophy conserving scheme robust to vanishing layers
      ! Note: Heffs are in lieu of h_at_u that should be returned by the
      !       continuity solver. AJA
      do j=jsq,jeq ; do i=is,ie
        heff1 = abs(uh(i,j,k) * g%IdyCu(i,j)) / (eps_vel+abs(u(i,j,k)))
        heff1 = max(heff1, min(h(i,j,k),h(i+1,j,k)))
        heff1 = min(heff1, max(h(i,j,k),h(i+1,j,k)))
        heff2 = abs(uh(i-1,j,k) * g%IdyCu(i-1,j)) / (eps_vel+abs(u(i-1,j,k)))
        heff2 = max(heff2, min(h(i-1,j,k),h(i,j,k)))
        heff2 = min(heff2, max(h(i-1,j,k),h(i,j,k)))
        heff3 = abs(uh(i,j+1,k) * g%IdyCu(i,j+1)) / (eps_vel+abs(u(i,j+1,k)))
        heff3 = max(heff3, min(h(i,j+1,k),h(i+1,j+1,k)))
        heff3 = min(heff3, max(h(i,j+1,k),h(i+1,j+1,k)))
        heff4 = abs(uh(i-1,j+1,k) * g%IdyCu(i-1,j+1)) / (eps_vel+abs(u(i-1,j+1,k)))
        heff4 = max(heff4, min(h(i-1,j+1,k),h(i,j+1,k)))
        heff4 = min(heff4, max(h(i-1,j+1,k),h(i,j+1,k)))
        if (cs%PV_Adv_Scheme == pv_adv_centered) then
          cav(i,j,k) = - 0.5*(abs_vort(i,j)+abs_vort(i-1,j)) * &
                         ((uh(i  ,j  ,k)+uh(i-1,j+1,k)) +      &
                          (uh(i-1,j  ,k)+uh(i  ,j+1,k)) ) /    &
                      (h_tiny + ((heff1+heff4) +(heff2+heff3)) ) * g%IdyCv(i,j)
        elseif (cs%PV_Adv_Scheme == pv_adv_upwind1) then
          uheff = ((uh(i  ,j  ,k)+uh(i-1,j+1,k)) +      &
                   (uh(i-1,j  ,k)+uh(i  ,j+1,k)) )
          quheff = 0.5*( (abs_vort(i,j)+abs_vort(i-1,j))*uheff &
                        -(abs_vort(i,j)-abs_vort(i-1,j))*abs(uheff) )
          cav(i,j,k) = - quheff / &
                       (h_tiny + ((heff1+heff4) +(heff2+heff3)) ) * g%IdyCv(i,j)
        endif
      enddo ; enddo
    endif
    ! Add in the additonal terms with Arakawa & Lamb.
    if ((cs%Coriolis_Scheme == arakawa_lamb81) .or. &
        (cs%Coriolis_Scheme == al_blend)) then ; do j=jsq,jeq ; do i=is,ie
      cav(i,j,k) = cav(i,j,k) + &
            (ep_v(i,j)*vh(i,j-1,k) - ep_v(i,j+1)*vh(i,j+1,k)) * g%IdyCv(i,j)
    enddo ; enddo ; endif

    if (cs%bound_Coriolis) then
      do j=jsq,jeq ; do i=is,ie
        max_fu = max(max_fuq(i,j),max_fuq(i-1,j))
        min_fu = min(min_fuq(i,j),min_fuq(i-1,j))
        if (cav(i,j,k) > max_fu) then
            cav(i,j,k) = max_fu
        else
          if (cav(i,j,k) < min_fu) cav(i,j,k) = min_fu
        endif
      enddo ; enddo
    endif

    ! Term - d(KE)/dy.
    do j=jsq,jeq ; do i=is,ie
      cav(i,j,k) = cav(i,j,k) - key(i,j)
      if (associated(ad%gradKEv)) ad%gradKEv(i,j,k) = -key(i,j)
    enddo ; enddo

    if (associated(ad%rv_x_u) .or. associated(ad%rv_x_v)) then
      ! Calculate the Coriolis-like acceleration due to relative vorticity.
      if (cs%Coriolis_Scheme == sadourny75_energy) then
        if (associated(ad%rv_x_u)) then
          do j=jsq,jeq ; do i=is,ie
            ad%rv_x_u(i,j,k) = - 0.25* &
              (q2(i-1,j)*(uh(i-1,j,k) + uh(i-1,j+1,k)) + &
               q2(i,j)*(uh(i,j,k) + uh(i,j+1,k))) * g%IdyCv(i,j)
          enddo ; enddo
        endif

        if (associated(ad%rv_x_v)) then
          do j=js,je ; do i=isq,ieq
            ad%rv_x_v(i,j,k) = 0.25 * &
              (q2(i,j) * (vh(i+1,j,k) + vh(i,j,k)) + &
               q2(i,j-1) * (vh(i,j-1,k) + vh(i+1,j-1,k))) * g%IdxCu(i,j)
          enddo ; enddo
        endif
      else
        if (associated(ad%rv_x_u)) then
          do j=jsq,jeq ; do i=is,ie
            ad%rv_x_u(i,j,k) = -g%IdyCv(i,j) * c1_12 * &
              ((q2(i,j) + q2(i-1,j) + q2(i-1,j-1)) * uh(i-1,j,k) + &
               (q2(i-1,j) + q2(i,j) + q2(i,j-1)) * uh(i,j,k) + &
               (q2(i-1,j) + q2(i,j+1) + q2(i,j)) * uh(i,j+1,k) + &
               (q2(i,j) + q2(i-1,j+1) + q2(i-1,j)) * uh(i-1,j+1,k))
          enddo ; enddo
        endif

        if (associated(ad%rv_x_v)) then
          do j=js,je ; do i=isq,ieq
            ad%rv_x_v(i,j,k) = g%IdxCu(i,j) * c1_12 * &
              ((q2(i+1,j) + q2(i,j) + q2(i,j-1)) * vh(i+1,j,k) + &
               (q2(i-1,j) + q2(i,j) + q2(i,j-1)) * vh(i,j,k) + &
               (q2(i-1,j-1) + q2(i,j) + q2(i,j-1)) * vh(i,j-1,k) + &
               (q2(i+1,j-1) + q2(i,j) + q2(i,j-1)) * vh(i+1,j-1,k))
          enddo ; enddo
        endif
      endif
    endif

  enddo ! k-loop.

  ! Here the various Coriolis-related derived quantities are offered for averaging.
  if (query_averaging_enabled(cs%diag)) then
    if (cs%id_rv > 0) call post_data(cs%id_rv, rv, cs%diag)
    if (cs%id_PV > 0) call post_data(cs%id_PV, pv, cs%diag)
    if (cs%id_gKEu>0) call post_data(cs%id_gKEu, ad%gradKEu, cs%diag)
    if (cs%id_gKEv>0) call post_data(cs%id_gKEv, ad%gradKEv, cs%diag)
    if (cs%id_rvxu > 0) call post_data(cs%id_rvxu, ad%rv_x_u, cs%diag)
    if (cs%id_rvxv > 0) call post_data(cs%id_rvxv, ad%rv_x_v, cs%diag)

    ! Diagnostics for terms multiplied by fractional thicknesses

    ! 3D diagnostics hf_gKEu etc. are commented because there is no clarity on proper remapping grid option.
    ! The code is retained for degugging purposes in the future.
    !if (CS%id_hf_gKEu > 0) then
    !  allocate(hf_gKEu(G%IsdB:G%IedB,G%jsd:G%jed,G%ke))
    !  do k=1,nz ; do j=js,je ; do I=Isq,Ieq
    !    hf_gKEu(I,j,k) = AD%gradKEu(I,j,k) * AD%diag_hfrac_u(I,j,k)
    !  enddo ; enddo ; enddo
    !  call post_data(CS%id_hf_gKEu, hf_gKEu, CS%diag)
    !endif

    !if (CS%id_hf_gKEv > 0) then
    !  allocate(hf_gKEv(G%isd:G%ied,G%JsdB:G%JedB,G%ke))
    !  do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
    !    hf_gKEv(i,J,k) = AD%gradKEv(i,J,k) * AD%diag_hfrac_v(i,J,k)
    !  enddo ; enddo ; enddo
    !  call post_data(CS%id_hf_gKEv, hf_gKEv, CS%diag)
    !endif

    if (cs%id_hf_gKEu_2d > 0) then
      allocate(hf_gkeu_2d(g%IsdB:g%IedB,g%jsd:g%jed))
      hf_gkeu_2d(:,:) = 0.0
      do k=1,nz ; do j=js,je ; do i=isq,ieq
        hf_gkeu_2d(i,j) = hf_gkeu_2d(i,j) + ad%gradKEu(i,j,k) * ad%diag_hfrac_u(i,j,k)
      enddo ; enddo ; enddo
      call post_data(cs%id_hf_gKEu_2d, hf_gkeu_2d, cs%diag)
      deallocate(hf_gkeu_2d)
    endif

    if (cs%id_hf_gKEv_2d > 0) then
      allocate(hf_gkev_2d(g%isd:g%ied,g%JsdB:g%JedB))
      hf_gkev_2d(:,:) = 0.0
      do k=1,nz ; do j=jsq,jeq ; do i=is,ie
        hf_gkev_2d(i,j) = hf_gkev_2d(i,j) + ad%gradKEv(i,j,k) * ad%diag_hfrac_v(i,j,k)
      enddo ; enddo ; enddo
      call post_data(cs%id_hf_gKEv_2d, hf_gkev_2d, cs%diag)
      deallocate(hf_gkev_2d)
    endif

    !if (CS%id_hf_rvxv > 0) then
    !  allocate(hf_rvxv(G%IsdB:G%IedB,G%jsd:G%jed,G%ke))
    !  do k=1,nz ; do j=js,je ; do I=Isq,Ieq
    !    hf_rvxv(I,j,k) = AD%rv_x_v(I,j,k) * AD%diag_hfrac_u(I,j,k)
    !  enddo ; enddo ; enddo
    !  call post_data(CS%id_hf_rvxv, hf_rvxv, CS%diag)
    !endif

    !if (CS%id_hf_rvxu > 0) then
    !  allocate(hf_rvxu(G%isd:G%ied,G%JsdB:G%JedB,G%ke))
    !  do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
    !    hf_rvxu(i,J,k) = AD%rv_x_u(i,J,k) * AD%diag_hfrac_v(i,J,k)
    !  enddo ; enddo ; enddo
    !  call post_data(CS%id_hf_rvxu, hf_rvxu, CS%diag)
    !endif

    if (cs%id_hf_rvxv_2d > 0) then
      allocate(hf_rvxv_2d(g%IsdB:g%IedB,g%jsd:g%jed))
      hf_rvxv_2d(:,:) = 0.0
      do k=1,nz ; do j=js,je ; do i=isq,ieq
        hf_rvxv_2d(i,j) = hf_rvxv_2d(i,j) + ad%rv_x_v(i,j,k) * ad%diag_hfrac_u(i,j,k)
      enddo ; enddo ; enddo
      call post_data(cs%id_hf_rvxv_2d, hf_rvxv_2d, cs%diag)
      deallocate(hf_rvxv_2d)
    endif

    if (cs%id_hf_rvxu_2d > 0) then
      allocate(hf_rvxu_2d(g%isd:g%ied,g%JsdB:g%JedB))
      hf_rvxu_2d(:,:) = 0.0
      do k=1,nz ; do j=jsq,jeq ; do i=is,ie
        hf_rvxu_2d(i,j) = hf_rvxu_2d(i,j) + ad%rv_x_u(i,j,k) * ad%diag_hfrac_v(i,j,k)
      enddo ; enddo ; enddo
      call post_data(cs%id_hf_rvxu_2d, hf_rvxu_2d, cs%diag)
      deallocate(hf_rvxu_2d)
    endif
  endif

end subroutine coradcalc


!> Calculates the acceleration due to the gradient of kinetic energy.
subroutine gradke(u, v, h, KE, KEx, KEy, k, OBC, G, US, CS)
  type(ocean_grid_type),                      intent(in)  :: G !< Ocen grid structure
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)),  intent(in)  :: u !< Zonal velocity [L T-1 ~> m s-1]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)),  intent(in)  :: v !< Meridional velocity [L T-1 ~> m s-1]
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),   intent(in)  :: h !< Layer thickness [H ~> m or kg m-2]
  real, dimension(SZI_(G) ,SZJ_(G) ),         intent(out) :: KE !< Kinetic energy per unit mass [L2 T-2 ~> m2 s-2]
  real, dimension(SZIB_(G),SZJ_(G) ),         intent(out) :: KEx !< Zonal acceleration due to kinetic
                                                                 !! energy gradient [L T-2 ~> m s-2]
  real, dimension(SZI_(G) ,SZJB_(G)),         intent(out) :: KEy !< Meridional acceleration due to kinetic
                                                                 !! energy gradient [L T-2 ~> m s-2]
  integer,                                    intent(in)  :: k !< Layer number to calculate for
  type(ocean_obc_type),                       pointer     :: OBC !< Open boundary control structure
  type(unit_scale_type),                      intent(in)  :: US  !< A dimensional unit scaling type
  type(coriolisadv_cs),                       pointer     :: CS !< Control structure for MOM_CoriolisAdv
  ! Local variables
  real :: um, up, vm, vp         ! Temporary variables [L T-1 ~> m s-1].
  real :: um2, up2, vm2, vp2     ! Temporary variables [L2 T-2 ~> m2 s-2].
  real :: um2a, up2a, vm2a, vp2a ! Temporary variables [L4 T-2 ~> m4 s-2].
  integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, n

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB


  ! Calculate KE (Kinetic energy for use in the -grad(KE) acceleration term).
  if (cs%KE_Scheme == ke_arakawa) then
    ! The following calculation of Kinetic energy includes the metric terms
    ! identified in Arakawa & Lamb 1982 as important for KE conservation.  It
    ! also includes the possibility of partially-blocked tracer cell faces.
    do j=jsq,jeq+1 ; do i=isq,ieq+1
      ke(i,j) = ( ( g%areaCu( i ,j)*(u( i ,j,k)*u( i ,j,k)) + &
                    g%areaCu(i-1,j)*(u(i-1,j,k)*u(i-1,j,k)) ) + &
                  ( g%areaCv(i, j )*(v(i, j ,k)*v(i, j ,k)) + &
                    g%areaCv(i,j-1)*(v(i,j-1,k)*v(i,j-1,k)) ) )*0.25*g%IareaT(i,j)
    enddo ; enddo
  elseif (cs%KE_Scheme == ke_simple_gudonov) then
    ! The following discretization of KE is based on the one-dimensinal Gudonov
    ! scheme which does not take into account any geometric factors
    do j=jsq,jeq+1 ; do i=isq,ieq+1
      up = 0.5*( u(i-1,j,k) + abs( u(i-1,j,k) ) ) ; up2 = up*up
      um = 0.5*( u( i ,j,k) - abs( u( i ,j,k) ) ) ; um2 = um*um
      vp = 0.5*( v(i,j-1,k) + abs( v(i,j-1,k) ) ) ; vp2 = vp*vp
      vm = 0.5*( v(i, j ,k) - abs( v(i, j ,k) ) ) ; vm2 = vm*vm
      ke(i,j) = ( max(up2,um2) + max(vp2,vm2) ) *0.5
    enddo ; enddo
  elseif (cs%KE_Scheme == ke_gudonov) then
    ! The following discretization of KE is based on the one-dimensinal Gudonov
    ! scheme but has been adapted to take horizontal grid factors into account
    do j=jsq,jeq+1 ; do i=isq,ieq+1
      up = 0.5*( u(i-1,j,k) + abs( u(i-1,j,k) ) ) ; up2a = up*up*g%areaCu(i-1,j)
      um = 0.5*( u( i ,j,k) - abs( u( i ,j,k) ) ) ; um2a = um*um*g%areaCu( i ,j)
      vp = 0.5*( v(i,j-1,k) + abs( v(i,j-1,k) ) ) ; vp2a = vp*vp*g%areaCv(i,j-1)
      vm = 0.5*( v(i, j ,k) - abs( v(i, j ,k) ) ) ; vm2a = vm*vm*g%areaCv(i, j )
      ke(i,j) = ( max(um2a,up2a) + max(vm2a,vp2a) )*0.5*g%IareaT(i,j)
    enddo ; enddo
  endif

  ! Term - d(KE)/dx.
  do j=js,je ; do i=isq,ieq
    kex(i,j) = (ke(i+1,j) - ke(i,j)) * g%IdxCu(i,j)
  enddo ; enddo

  ! Term - d(KE)/dy.
  do j=jsq,jeq ; do i=is,ie
    key(i,j) = (ke(i,j+1) - ke(i,j)) * g%IdyCv(i,j)
  enddo ; enddo

  if (associated(obc)) then
    do n=1,obc%number_of_segments
      if (obc%segment(n)%is_N_or_S) then
        do i=obc%segment(n)%HI%isd,obc%segment(n)%HI%ied
          key(i,obc%segment(n)%HI%JsdB) = 0.
        enddo
      elseif (obc%segment(n)%is_E_or_W) then
        do j=obc%segment(n)%HI%jsd,obc%segment(n)%HI%jed
          kex(obc%segment(n)%HI%IsdB,j) = 0.
        enddo
      endif
    enddo
  endif

end subroutine gradke

!> Initializes the control structure for coriolisadv_cs
subroutine coriolisadv_init(Time, G, GV, US, param_file, diag, AD, CS)
  type(time_type), target, intent(in)    :: time !< Current model time
  type(ocean_grid_type),   intent(in)    :: g  !< Ocean grid structure
  type(verticalgrid_type), intent(in)    :: gv !< Vertical grid structure
  type(unit_scale_type),   intent(in)    :: us  !< A dimensional unit scaling type
  type(param_file_type),   intent(in)    :: param_file !< Runtime parameter handles
  type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure
  type(accel_diag_ptrs),   target, intent(inout) :: ad !< Strorage for acceleration diagnostics
  type(coriolisadv_cs),    pointer       :: cs !< Control structure fro MOM_CoriolisAdv
  ! Local variables
! This include declares and sets the variable "version".
#include "version_variable.h"
  character(len=40)  :: mdl = "MOM_CoriolisAdv" ! This module's name.
  character(len=20)  :: tmpstr
  character(len=400) :: mesg
  integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb, nz

  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = gv%ke
  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB

  if (associated(cs)) then
    call mom_error(warning, "CoriolisAdv_init called with associated control structure.")
    return
  endif
  allocate(cs)

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

  ! Read all relevant parameters and write them to the model log.
  call log_version(param_file, mdl, version, "")
  call get_param(param_file, mdl, "NOSLIP", cs%no_slip, &
                 "If true, no slip boundary conditions are used; otherwise "//&
                 "free slip boundary conditions are assumed. The "//&
                 "implementation of the free slip BCs on a C-grid is much "//&
                 "cleaner than the no slip BCs. The use of free slip BCs "//&
                 "is strongly encouraged, and no slip BCs are not used with "//&
                 "the biharmonic viscosity.", default=.false.)

  call get_param(param_file, mdl, "CORIOLIS_EN_DIS", cs%Coriolis_En_Dis, &
                 "If true, two estimates of the thickness fluxes are used "//&
                 "to estimate the Coriolis term, and the one that "//&
                 "dissipates energy relative to the other one is used.", &
                 default=.false.)

  ! Set %Coriolis_Scheme
  ! (Select the baseline discretization for the Coriolis term)
  call get_param(param_file, mdl, "CORIOLIS_SCHEME", tmpstr, &
                 "CORIOLIS_SCHEME selects the discretization for the "//&
                 "Coriolis terms. Valid values are: \n"//&
                 "\t SADOURNY75_ENERGY - Sadourny, 1975; energy cons. \n"//&
                 "\t ARAKAWA_HSU90     - Arakawa & Hsu, 1990 \n"//&
                 "\t SADOURNY75_ENSTRO - Sadourny, 1975; enstrophy cons. \n"//&
                 "\t ARAKAWA_LAMB81    - Arakawa & Lamb, 1981; En. + Enst.\n"//&
                 "\t ARAKAWA_LAMB_BLEND - A blend of Arakawa & Lamb with \n"//&
                 "\t                      Arakawa & Hsu and Sadourny energy", &
                 default=sadourny75_energy_string)
  tmpstr = uppercase(tmpstr)
  select case (tmpstr)
    case (sadourny75_energy_string)
      cs%Coriolis_Scheme = sadourny75_energy
    case (arakawa_hsu_string)
      cs%Coriolis_Scheme = arakawa_hsu90
    case (sadourny75_enstro_string)
      cs%Coriolis_Scheme = sadourny75_enstro
    case (arakawa_lamb_string)
      cs%Coriolis_Scheme = arakawa_lamb81
    case (al_blend_string)
      cs%Coriolis_Scheme = al_blend
    case (robust_enstro_string)
      cs%Coriolis_Scheme = robust_enstro
      cs%Coriolis_En_Dis = .false.
    case default
      call mom_mesg('CoriolisAdv_init: Coriolis_Scheme ="'//trim(tmpstr)//'"', 0)
      call mom_error(fatal, "CoriolisAdv_init: Unrecognized setting "// &
            "#define CORIOLIS_SCHEME "//trim(tmpstr)//" found in input file.")
  end select
  if (cs%Coriolis_Scheme == al_blend) then
    call get_param(param_file, mdl, "CORIOLIS_BLEND_WT_LIN", cs%wt_lin_blend, &
                 "A weighting value for the ratio of inverse thicknesses, "//&
                 "beyond which the blending between Sadourny Energy and "//&
                 "Arakawa & Hsu goes linearly to 0 when CORIOLIS_SCHEME "//&
                 "is ARAWAKA_LAMB_BLEND. This must be between 1 and 1e-16.", &
                 units="nondim", default=0.125)
    call get_param(param_file, mdl, "CORIOLIS_BLEND_F_EFF_MAX", cs%F_eff_max_blend, &
                 "The factor by which the maximum effective Coriolis "//&
                 "acceleration from any point can be increased when "//&
                 "blending different discretizations with the "//&
                 "ARAKAWA_LAMB_BLEND Coriolis scheme.  This must be "//&
                 "greater than 2.0 (the max value for Sadourny energy).", &
                 units="nondim", default=4.0)
    cs%wt_lin_blend = min(1.0, max(cs%wt_lin_blend,1e-16))
    if (cs%F_eff_max_blend < 2.0) call mom_error(warning, "CoriolisAdv_init: "//&
           "CORIOLIS_BLEND_F_EFF_MAX should be at least 2.")
  endif

  mesg = "If true, the Coriolis terms at u-points are bounded by "//&
         "the four estimates of (f+rv)v from the four neighboring "//&
         "v-points, and similarly at v-points."
  if (cs%Coriolis_En_Dis .and. (cs%Coriolis_Scheme == sadourny75_energy)) then
    mesg = trim(mesg)//"  This option is "//&
                 "always effectively false with CORIOLIS_EN_DIS defined and "//&
                 "CORIOLIS_SCHEME set to "//trim(sadourny75_energy_string)//"."
  else
    mesg = trim(mesg)//"  This option would "//&
                 "have no effect on the SADOURNY Coriolis scheme if it "//&
                 "were possible to use centered difference thickness fluxes."
  endif
  call get_param(param_file, mdl, "BOUND_CORIOLIS", cs%bound_Coriolis, mesg, &
                 default=.false.)
  if ((cs%Coriolis_En_Dis .and. (cs%Coriolis_Scheme == sadourny75_energy)) .or. &
      (cs%Coriolis_Scheme == robust_enstro)) cs%bound_Coriolis = .false.

  ! Set KE_Scheme (selects discretization of KE)
  call get_param(param_file, mdl, "KE_SCHEME", tmpstr, &
                 "KE_SCHEME selects the discretization for acceleration "//&
                 "due to the kinetic energy gradient. Valid values are: \n"//&
                 "\t KE_ARAKAWA, KE_SIMPLE_GUDONOV, KE_GUDONOV", &
                 default=ke_arakawa_string)
  tmpstr = uppercase(tmpstr)
  select case (tmpstr)
    case (ke_arakawa_string); cs%KE_Scheme = ke_arakawa
    case (ke_simple_gudonov_string); cs%KE_Scheme = ke_simple_gudonov
    case (ke_gudonov_string); cs%KE_Scheme = ke_gudonov
    case default
      call mom_mesg('CoriolisAdv_init: KE_Scheme ="'//trim(tmpstr)//'"', 0)
      call mom_error(fatal, "CoriolisAdv_init: "// &
               "#define KE_SCHEME "//trim(tmpstr)//" in input file is invalid.")
  end select

  ! Set PV_Adv_Scheme (selects discretization of PV advection)
  call get_param(param_file, mdl, "PV_ADV_SCHEME", tmpstr, &
                 "PV_ADV_SCHEME selects the discretization for PV "//&
                 "advection. Valid values are: \n"//&
                 "\t PV_ADV_CENTERED - centered (aka Sadourny, 75) \n"//&
                 "\t PV_ADV_UPWIND1  - upwind, first order", &
                 default=pv_adv_centered_string)
  select case (uppercase(tmpstr))
    case (pv_adv_centered_string)
      cs%PV_Adv_Scheme = pv_adv_centered
    case (pv_adv_upwind1_string)
      cs%PV_Adv_Scheme = pv_adv_upwind1
    case default
      call mom_mesg('CoriolisAdv_init: PV_Adv_Scheme ="'//trim(tmpstr)//'"', 0)
      call mom_error(fatal, "CoriolisAdv_init: "// &
                     "#DEFINE PV_ADV_SCHEME in input file is invalid.")
  end select

  cs%id_rv = register_diag_field('ocean_model', 'RV', diag%axesBL, time, &
     'Relative Vorticity', 's-1', conversion=us%s_to_T)

  cs%id_PV = register_diag_field('ocean_model', 'PV', diag%axesBL, time, &
     'Potential Vorticity', 'm-1 s-1', conversion=gv%m_to_H*us%s_to_T)

  cs%id_gKEu = register_diag_field('ocean_model', 'gKEu', diag%axesCuL, time, &
     'Zonal Acceleration from Grad. Kinetic Energy', 'm-1 s-2', conversion=us%L_T2_to_m_s2)
  if (cs%id_gKEu > 0) call safe_alloc_ptr(ad%gradKEu,isdb,iedb,jsd,jed,nz)

  cs%id_gKEv = register_diag_field('ocean_model', 'gKEv', diag%axesCvL, time, &
     'Meridional Acceleration from Grad. Kinetic Energy', 'm-1 s-2', conversion=us%L_T2_to_m_s2)
  if (cs%id_gKEv > 0) call safe_alloc_ptr(ad%gradKEv,isd,ied,jsdb,jedb,nz)

  cs%id_rvxu = register_diag_field('ocean_model', 'rvxu', diag%axesCvL, time, &
     'Meridional Acceleration from Relative Vorticity', 'm-1 s-2', conversion=us%L_T2_to_m_s2)
  if (cs%id_rvxu > 0) call safe_alloc_ptr(ad%rv_x_u,isd,ied,jsdb,jedb,nz)

  cs%id_rvxv = register_diag_field('ocean_model', 'rvxv', diag%axesCuL, time, &
     'Zonal Acceleration from Relative Vorticity', 'm-1 s-2', conversion=us%L_T2_to_m_s2)
  if (cs%id_rvxv > 0) call safe_alloc_ptr(ad%rv_x_v,isdb,iedb,jsd,jed,nz)

  !CS%id_hf_gKEu = register_diag_field('ocean_model', 'hf_gKEu', diag%axesCuL, Time, &
  !   'Fractional Thickness-weighted Zonal Acceleration from Grad. Kinetic Energy', &
  !   'm-1 s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)
  !if (CS%id_hf_gKEu > 0) then
  !  call safe_alloc_ptr(AD%gradKEu,IsdB,IedB,jsd,jed,nz)
  !  call safe_alloc_ptr(AD%diag_hfrac_u,IsdB,IedB,jsd,jed,nz)
  !endif

  !CS%id_hf_gKEv = register_diag_field('ocean_model', 'hf_gKEv', diag%axesCvL, Time, &
  !   'Fractional Thickness-weighted Meridional Acceleration from Grad. Kinetic Energy', &
  !   'm-1 s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)
  !if (CS%id_hf_gKEv > 0) then
  !  call safe_alloc_ptr(AD%gradKEv,isd,ied,JsdB,JedB,nz)
  !  call safe_alloc_ptr(AD%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
  !endif

  cs%id_hf_gKEu_2d = register_diag_field('ocean_model', 'hf_gKEu_2d', diag%axesCu1, time, &
     'Depth-sum Fractional Thickness-weighted Zonal Acceleration from Grad. Kinetic Energy', &
     'm-1 s-2', conversion=us%L_T2_to_m_s2)
  if (cs%id_hf_gKEu_2d > 0) then
    call safe_alloc_ptr(ad%gradKEu,isdb,iedb,jsd,jed,nz)
    call safe_alloc_ptr(ad%diag_hfrac_u,isdb,iedb,jsd,jed,nz)
  endif

  cs%id_hf_gKEv_2d = register_diag_field('ocean_model', 'hf_gKEv_2d', diag%axesCv1, time, &
     'Depth-sum Fractional Thickness-weighted Meridional Acceleration from Grad. Kinetic Energy', &
     'm-1 s-2', conversion=us%L_T2_to_m_s2)
  if (cs%id_hf_gKEv_2d > 0) then
    call safe_alloc_ptr(ad%gradKEv,isd,ied,jsdb,jedb,nz)
    call safe_alloc_ptr(ad%diag_hfrac_v,isd,ied,jsd,jedb,nz)
  endif

  !CS%id_hf_rvxu = register_diag_field('ocean_model', 'hf_rvxu', diag%axesCvL, Time, &
  !   'Fractional Thickness-weighted Meridional Acceleration from Relative Vorticity', &
  !   'm-1 s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)
  !if (CS%id_hf_rvxu > 0) then
  !  call safe_alloc_ptr(AD%rv_x_u,isd,ied,JsdB,JedB,nz)
  !  call safe_alloc_ptr(AD%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
  !endif

  !CS%id_hf_rvxv = register_diag_field('ocean_model', 'hf_rvxv', diag%axesCuL, Time, &
  !   'Fractional Thickness-weighted Zonal Acceleration from Relative Vorticity', &
  !   'm-1 s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)
  !if (CS%id_hf_rvxv > 0) then
  !  call safe_alloc_ptr(AD%rv_x_v,IsdB,IedB,jsd,jed,nz)
  !  call safe_alloc_ptr(AD%diag_hfrac_u,IsdB,IedB,jsd,jed,nz)
  !endif

  cs%id_hf_rvxu_2d = register_diag_field('ocean_model', 'hf_rvxu_2d', diag%axesCv1, time, &
     'Fractional Thickness-weighted Meridional Acceleration from Relative Vorticity', &
     'm-1 s-2', conversion=us%L_T2_to_m_s2)
  if (cs%id_hf_rvxu_2d > 0) then
    call safe_alloc_ptr(ad%rv_x_u,isd,ied,jsdb,jedb,nz)
    call safe_alloc_ptr(ad%diag_hfrac_v,isd,ied,jsd,jedb,nz)
  endif

  cs%id_hf_rvxv_2d = register_diag_field('ocean_model', 'hf_rvxv_2d', diag%axesCu1, time, &
     'Depth-sum Fractional Thickness-weighted Zonal Acceleration from Relative Vorticity', &
     'm-1 s-2', conversion=us%L_T2_to_m_s2)
  if (cs%id_hf_rvxv_2d > 0) then
    call safe_alloc_ptr(ad%rv_x_v,isdb,iedb,jsd,jed,nz)
    call safe_alloc_ptr(ad%diag_hfrac_u,isdb,iedb,jsd,jed,nz)
  endif

end subroutine coriolisadv_init

!> Destructor for coriolisadv_cs
subroutine coriolisadv_end(CS)
  type(coriolisadv_cs), pointer :: cs !< Control structure fro MOM_CoriolisAdv
  deallocate(cs)
end subroutine coriolisadv_end

!> \namespace mom_coriolisadv
!!
!! This file contains the subroutine that calculates the time
!! derivatives of the velocities due to Coriolis acceleration and
!! momentum advection.  This subroutine uses either a vorticity
!! advection scheme from Arakawa and Hsu, Mon. Wea. Rev. 1990, or
!! Sadourny's (JAS 1975) energy conserving scheme.  Both have been
!! modified to use general orthogonal coordinates as described in
!! Arakawa and Lamb, Mon. Wea. Rev. 1981.  Both schemes are second
!! order accurate, and allow for vanishingly small layer thicknesses.
!! The Arakawa and Hsu scheme globally conserves both total energy
!! and potential enstrophy in the limit of nondivergent flow.
!! Sadourny's energy conserving scheme conserves energy if the flow
!! is nondivergent or centered difference thickness fluxes are used.
!!
!! A small fragment of the grid is shown below:
!! \verbatim
!!
!!    j+1  x ^ x ^ x   At x:  q, CoriolisBu
!!    j+1  > o > o >   At ^:  v, CAv, vh
!!    j    x ^ x ^ x   At >:  u, CAu, uh, a, b, c, d
!!    j    > o > o >   At o:  h, KE
!!    j-1  x ^ x ^ x
!!        i-1  i  i+1  At x & ^:
!!           i  i+1    At > & o:
!! \endverbatim
!!
!! The boundaries always run through q grid points (x).

end module mom_coriolisadv