MOM_dynamics_unsplit_RK2.F90

!> Time steps the ocean dynamics with an unsplit quasi 2nd order Runge-Kutta scheme
module mom_dynamics_unsplit_rk2

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

!********+*********+*********+*********+*********+*********+*********+**
!*                                                                     *
!*  By Alistair Adcroft and Robert Hallberg, 2010-2012                 *
!*                                                                     *
!*    This file contains code that does the time-stepping of the       *
!*  adiabatic dynamic core, in this case with a pseudo-second order    *
!*  Runge-Kutta time stepping scheme for the momentum and a forward-   *
!*  backward coupling between the momentum and continuity equations,   *
!*  but without any splitting between the baroclinic and barotropic    *
!*  modes. Apart from the lack of splitting, this is closely analogous *
!*  to the split time stepping scheme, and efforts have been taken to  *
!*  ensure that for certain configurations (e.g., very short           *
!*  baroclinic time steps, a single barotropic step per baroclinic     *
!*  step, and particular choices about how to coupled the baroclinic   *
!*  and barotropic solves, the two solutions reproduce each other.     *
!*  Although this time stepping scheme is not very efficient with a    *
!*  large number of layers, it is valuable for verifying the proper    *
!*  behavior of the more complicated split time stepping scheme, and   *
!*  is not too inefficient for use with only a few layers.             *
!*                                                                     *
!*    The subroutine step_MOM_dyn_unsplit_RK2 actually does the time   *
!*  stepping, while register_restarts_dyn_unsplit_RK2 sets the fields  *
!*  that are found in a full restart file with this scheme, and        *
!*  initialize_dyn_unsplit_RK2 initializes the cpu clocks that are     *
!*  used in this module.  For largely historical reasons, this module  *
!*  does not have its own control structure, but shares the same       *
!*  control structure with MOM.F90 and the other MOM_dynamics_...      *
!*  modules.                                                           *
!*                                                                     *
!*  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, tauy, vbt, vhtr  *
!*    j    x ^ x ^ x   At >:  u, PFu, CAu, uh, diffu, taux, ubt, uhtr  *
!*    j    > o > o >   At o:  h, bathyT, eta, T, S, tr                 *
!*    j-1  x ^ x ^ x                                                   *
!*        i-1  i  i+1                                                  *
!*           i  i+1                                                    *
!*                                                                     *
!*  The boundaries always run through q grid points (x).               *
!*                                                                     *
!********+*********+*********+*********+*********+*********+*********+**

use mom_variables, only : vertvisc_type, thermo_var_ptrs
use mom_variables, only : ocean_internal_state, accel_diag_ptrs, cont_diag_ptrs
use mom_forcing_type, only : mech_forcing
use mom_checksum_packages, only : mom_thermo_chksum, mom_state_chksum, mom_accel_chksum
use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
use mom_cpu_clock, only : clock_component, clock_subcomponent
use mom_cpu_clock, only : clock_module_driver, clock_module, clock_routine
use mom_diag_mediator, only : diag_mediator_init, enable_averages
use mom_diag_mediator, only : disable_averaging, post_data, safe_alloc_ptr
use mom_diag_mediator, only : register_diag_field, register_static_field
use mom_diag_mediator, only : set_diag_mediator_grid, diag_ctrl
use mom_domains, only : mom_domains_init, pass_var, pass_vector
use mom_domains, only : pass_var_start, pass_var_complete
use mom_domains, only : pass_vector_start, pass_vector_complete
use mom_domains, only : to_south, to_west, to_all, cgrid_ne, scalar_pair
use mom_error_handler, only : mom_error, mom_mesg, fatal, warning, is_root_pe
use mom_error_handler, only : mom_set_verbosity
use mom_file_parser, only : get_param, log_version, param_file_type
use mom_get_input, only : directories
use mom_io, only : mom_io_init
use mom_restart, only : register_restart_field, query_initialized, save_restart
use mom_restart, only : restart_init, mom_restart_cs
use mom_time_manager, only : time_type, time_type_to_real, operator(+)
use mom_time_manager, only : operator(-), operator(>), operator(*), operator(/)

use mom_ale, only : ale_cs
use mom_boundary_update, only : update_obc_data, update_obc_cs
use mom_barotropic, only : barotropic_cs
use mom_continuity, only : continuity, continuity_init, continuity_cs, continuity_stencil
use mom_coriolisadv, only : coradcalc, coriolisadv_init, coriolisadv_cs
use mom_debugging, only : check_redundant
use mom_grid, only : ocean_grid_type
use mom_hor_index, only : hor_index_type
use mom_hor_visc, only : horizontal_viscosity, hor_visc_init, hor_visc_cs
use mom_lateral_mixing_coeffs, only : varmix_cs
use mom_meke_types, only : meke_type
use mom_open_boundary, only : ocean_obc_type
use mom_open_boundary, only : radiation_open_bdry_conds
use mom_open_boundary, only : open_boundary_zero_normal_flow
use mom_pressureforce, only : pressureforce, pressureforce_init, pressureforce_cs
use mom_set_visc, only : set_viscous_ml, set_visc_cs
use mom_thickness_diffuse, only : thickness_diffuse_cs
use mom_tidal_forcing, only : tidal_forcing_init, tidal_forcing_cs
use mom_unit_scaling, only : unit_scale_type
use mom_vert_friction, only : vertvisc, vertvisc_coef, vertvisc_init, vertvisc_cs
use mom_verticalgrid, only : verticalgrid_type, get_thickness_units
use mom_verticalgrid, only : get_flux_units, get_tr_flux_units

implicit none ; private

#include <MOM_memory.h>

!> MOM_dynamics_unsplit_RK2 module control structure
type, public :: mom_dyn_unsplit_rk2_cs ; private
  real allocable_, dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: &
    cau, &    !< CAu = f*v - u.grad(u) [L T-2 ~> m s-2].
    pfu, &    !< PFu = -dM/dx [L T-2 ~> m s-2].
    diffu     !< Zonal acceleration due to convergence of the along-isopycnal stress tensor [L T-2 ~> m s-2].

  real allocable_, dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: &
    cav, &    !< CAv = -f*u - u.grad(v) [L T-2 ~> m s-2].
    pfv, &    !< PFv = -dM/dy [L T-2 ~> m s-2].
    diffv     !< Meridional acceleration due to convergence of the along-isopycnal stress tensor [L T-2 ~> m s-2].

  real, pointer, dimension(:,:) :: taux_bot => null() !< frictional x-bottom stress from the ocean
                                                      !! to the seafloor [R L Z T-2 ~> Pa]
  real, pointer, dimension(:,:) :: tauy_bot => null() !< frictional y-bottom stress from the ocean
                                                      !! to the seafloor [R L Z T-2 ~> Pa]

  real    :: be      !< A nondimensional number from 0.5 to 1 that controls
                     !! the backward weighting of the time stepping scheme.
  real    :: begw    !< A nondimensional number from 0 to 1 that controls
                     !! the extent to which the treatment of gravity waves
                     !! is forward-backward (0) or simulated backward
                     !! Euler (1).  0 is almost always used.
  logical :: use_correct_dt_visc !< If true, use the correct timestep in the calculation of the
                                 !! turbulent mixed layer properties for viscosity.
                                 !! The default should be true, but it is false.
  logical :: debug   !< If true, write verbose checksums for debugging purposes.

  logical :: module_is_initialized = .false. !< Record whether this mouled has been initialzed.

  !>@{ Diagnostic IDs
  integer :: id_uh = -1, id_vh = -1
  integer :: id_pfu = -1, id_pfv = -1, id_cau = -1, id_cav = -1
  !>@}

  type(diag_ctrl), pointer :: diag => null() !< A structure that is used to
                                   !! regulate the timing of diagnostic output.
  type(accel_diag_ptrs), pointer :: adp => null() !< A structure pointing to the
                                   !! accelerations in the momentum equations,
                                   !! which can later be used to calculate
                                   !! derived diagnostics like energy budgets.
  type(cont_diag_ptrs), pointer :: cdp => null() !< A structure with pointers to
                                   !! various terms in the continuity equations,
                                   !! which can later be used to calculate
                                   !! derived diagnostics like energy budgets.

  ! The remainder of the structure points to child subroutines' control structures.
  !> A pointer to the horizontal viscosity control structure
  type(hor_visc_cs), pointer :: hor_visc_csp => null()
  !> A pointer to the continuity control structure
  type(continuity_cs), pointer :: continuity_csp => null()
  !> A pointer to the CoriolisAdv control structure
  type(coriolisadv_cs), pointer :: coriolisadv_csp => null()
  !> A pointer to the PressureForce control structure
  type(pressureforce_cs), pointer :: pressureforce_csp => null()
  !> A pointer to the vertvisc control structure
  type(vertvisc_cs), pointer :: vertvisc_csp => null()
  !> A pointer to the set_visc control structure
  type(set_visc_cs), pointer :: set_visc_csp => null()
  !> A pointer to the tidal forcing control structure
  type(tidal_forcing_cs), pointer :: tides_csp => null()
  !> A pointer to the ALE control structure.
  type(ale_cs), pointer :: ale_csp => null()

  type(ocean_obc_type), pointer :: obc => null() !< A pointer to an open boundary
     !! condition type that specifies whether, where, and what open boundary
     !! conditions are used.  If no open BCs are used, this pointer stays
     !! nullified.  Flather OBCs use open boundary_CS as well.
  !> A pointer to the update_OBC control structure
  type(update_obc_cs), pointer :: update_obc_csp => null()

end type mom_dyn_unsplit_rk2_cs


public step_mom_dyn_unsplit_rk2, register_restarts_dyn_unsplit_rk2
public initialize_dyn_unsplit_rk2, end_dyn_unsplit_rk2

!>@{ CPU time clock IDs
integer :: id_clock_cor, id_clock_pres, id_clock_vertvisc
integer :: id_clock_horvisc, id_clock_continuity, id_clock_mom_update
integer :: id_clock_pass, id_clock_pass_init
!>@}

contains

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

!> Step the MOM6 dynamics using an unsplit quasi-2nd order Runge-Kutta scheme
subroutine step_mom_dyn_unsplit_rk2(u_in, v_in, h_in, tv, visc, Time_local, dt, forces, &
                  p_surf_begin, p_surf_end, uh, vh, uhtr, vhtr, eta_av, G, GV, US, CS, &
                  VarMix, MEKE)
  type(ocean_grid_type),             intent(inout) :: G       !< The ocean's grid structure.
  type(verticalgrid_type),           intent(in)    :: GV      !< The ocean's vertical grid structure.
  type(unit_scale_type),             intent(in)    :: US      !< A dimensional unit scaling type
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: u_in !< The input and output zonal
                                                              !! velocity [L T-1 ~> m s-1].
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: v_in !< The input and output meridional
                                                              !! velocity [L T-1 ~> m s-1].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(inout) :: h_in !< The input and output layer thicknesses,
                                                              !! [H ~> m or kg m-2], depending on whether
                                                              !! the Boussinesq approximation is made.
  type(thermo_var_ptrs),             intent(in)    :: tv      !< A structure pointing to various
                                                              !! thermodynamic variables.
  type(vertvisc_type),               intent(inout) :: visc    !< A structure containing vertical
                                                              !! viscosities, bottom drag
                                                              !! viscosities, and related fields.
  type(time_type),                   intent(in)    :: Time_local   !< The model time at the end of
                                                              !! the time step.
  real,                              intent(in)    :: dt      !< The baroclinic dynamics time step [T ~> s].
  type(mech_forcing),                intent(in)    :: forces  !< A structure with the driving mechanical forces
  real, dimension(:,:),              pointer       :: p_surf_begin !< A pointer (perhaps NULL) to
                                                              !! the surface pressure at the beginning
                                                              !! of this dynamic step [R L2 T-2 ~> Pa].
  real, dimension(:,:),              pointer       :: p_surf_end   !< A pointer (perhaps NULL) to
                                                              !! the surface pressure at the end of
                                                              !! this dynamic step [R L2 T-2 ~> Pa].
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: uh !< The zonal volume or mass transport
                                                              !! [H L2 T-1 ~> m3 s-1 or kg s-1].
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: vh  !< The meridional volume or mass
                                                              !! transport [H L2 T-1 ~> m3 s-1 or kg s-1].
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: uhtr !< The accumulated zonal volume or
                                                              !! mass transport since the last
                                                              !! tracer advection [H L2 ~> m3 or kg].
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: vhtr !< The accumulated meridional volume
                                                              !! or mass transport since the last
                                                              !! tracer advection [H L2 ~> m3 or kg].
  real, dimension(SZI_(G),SZJ_(G)),  intent(out)   :: eta_av  !< The time-mean free surface height
                                                              !! or column mass [H ~> m or kg m-2].
  type(mom_dyn_unsplit_rk2_cs),      pointer       :: CS      !< The control structure set up by
                                                              !! initialize_dyn_unsplit_RK2.
  type(varmix_cs),                   pointer       :: VarMix  !< A pointer to a structure with
                                                              !! fields that specify the spatially
                                                              !! variable viscosities.
  type(meke_type),                   pointer       :: MEKE    !< A pointer to a structure containing
                                                              !! fields related to the Mesoscale
                                                              !! Eddy Kinetic Energy.
  ! Local variables
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_av, hp
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: up ! Predicted zonal velocities [L T-1 ~> m s-1]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: vp ! Predicted meridional velocities [L T-1 ~> m s-1]
  real, dimension(:,:), pointer :: p_surf => null()
  real :: dt_pred   ! The time step for the predictor part of the baroclinic time stepping [T ~> s]
  real :: dt_visc   ! The time step for a part of the update due to viscosity [T ~> s]
  logical :: dyn_p_surf
  integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
  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
  dt_pred = dt * cs%BE

  h_av(:,:,:) = 0; hp(:,:,:) = 0
  up(:,:,:) = 0
  vp(:,:,:) = 0

  dyn_p_surf = associated(p_surf_begin) .and. associated(p_surf_end)
  if (dyn_p_surf) then
    call safe_alloc_ptr(p_surf,g%isd,g%ied,g%jsd,g%jed) ; p_surf(:,:) = 0.0
  else
    p_surf => forces%p_surf
  endif

! Runge-Kutta second order accurate two step scheme is used to step
! all of the fields except h.  h is stepped separately.

  if (cs%debug) then
    call mom_state_chksum("Start Predictor ", u_in, v_in, h_in, uh, vh, g, gv, us)
  endif

! diffu = horizontal viscosity terms (u,h)
  call enable_averages(dt,time_local, cs%diag)
  call cpu_clock_begin(id_clock_horvisc)
  call horizontal_viscosity(u_in, v_in, h_in, cs%diffu, cs%diffv, meke, varmix, &
                            g, gv, us, cs%hor_visc_CSp)
  call cpu_clock_end(id_clock_horvisc)
  call disable_averaging(cs%diag)
  call pass_vector(cs%diffu, cs%diffv, g%Domain, clock=id_clock_pass)

! This continuity step is solely for the Coroilis terms, specifically in the
! denominator of PV and in the mass transport or PV.
! uh = u[n-1]*h[n-1/2]
! hp = h[n-1/2] + dt/2 div . uh
  call cpu_clock_begin(id_clock_continuity)
  ! This is a duplicate calculation of the last continuity from the previous step
  ! and could/should be optimized out. -AJA
  call continuity(u_in, v_in, h_in, hp, uh, vh, dt_pred, g, gv, us, cs%continuity_CSp, obc=cs%OBC)
  call cpu_clock_end(id_clock_continuity)
  call pass_var(hp, g%Domain, clock=id_clock_pass)
  call pass_vector(uh, vh, g%Domain, clock=id_clock_pass)

! h_av = (h + hp)/2  (used in PV denominator)
  call cpu_clock_begin(id_clock_mom_update)
  do k=1,nz ; do j=js-2,je+2 ; do i=is-2,ie+2
    h_av(i,j,k) = (h_in(i,j,k) + hp(i,j,k)) * 0.5
  enddo ; enddo ; enddo
  call cpu_clock_end(id_clock_mom_update)

! CAu = -(f+zeta)/h_av vh + d/dx KE  (function of u[n-1] and uh[n-1])
  call cpu_clock_begin(id_clock_cor)
  call coradcalc(u_in, v_in, h_av, uh, vh, cs%CAu, cs%CAv, cs%OBC, cs%ADp, &
                 g, gv, us, cs%CoriolisAdv_CSp)
  call cpu_clock_end(id_clock_cor)

! PFu = d/dx M(h_av,T,S)  (function of h[n-1/2])
  call cpu_clock_begin(id_clock_pres)
  if (dyn_p_surf) then ; do j=js-2,je+2 ; do i=is-2,ie+2
    p_surf(i,j) = 0.5*p_surf_begin(i,j) + 0.5*p_surf_end(i,j)
  enddo ; enddo ; endif
  call pressureforce(h_in, tv, cs%PFu, cs%PFv, g, gv, us, cs%PressureForce_CSp, cs%ALE_CSp, p_surf)
  call cpu_clock_end(id_clock_pres)
  call pass_vector(cs%PFu, cs%PFv, g%Domain, clock=id_clock_pass)
  call pass_vector(cs%CAu, cs%CAv, g%Domain, clock=id_clock_pass)

  if (associated(cs%OBC)) then; if (cs%OBC%update_OBC) then
    call update_obc_data(cs%OBC, g, gv, us, tv, h_in, cs%update_OBC_CSp, time_local)
  endif; endif
  if (associated(cs%OBC)) then
    call open_boundary_zero_normal_flow(cs%OBC, g, cs%PFu, cs%PFv)
    call open_boundary_zero_normal_flow(cs%OBC, g, cs%CAu, cs%CAv)
    call open_boundary_zero_normal_flow(cs%OBC, g, cs%diffu, cs%diffv)
  endif

! up+[n-1/2] = u[n-1] + dt_pred * (PFu + CAu)
  call cpu_clock_begin(id_clock_mom_update)
  do k=1,nz ; do j=js,je ; do i=isq,ieq
    up(i,j,k) = g%mask2dCu(i,j) * (u_in(i,j,k) + dt_pred * &
                   ((cs%PFu(i,j,k) + cs%CAu(i,j,k)) + cs%diffu(i,j,k)))
  enddo ; enddo ; enddo
  do k=1,nz ; do j=jsq,jeq ; do i=is,ie
    vp(i,j,k) = g%mask2dCv(i,j) * (v_in(i,j,k) + dt_pred * &
                   ((cs%PFv(i,j,k) + cs%CAv(i,j,k)) + cs%diffv(i,j,k)))
  enddo ; enddo ; enddo
  call cpu_clock_end(id_clock_mom_update)

  if (cs%debug) &
    call mom_accel_chksum("Predictor 1 accel", cs%CAu, cs%CAv, cs%PFu, cs%PFv,&
                          cs%diffu, cs%diffv, g, gv, us)

 ! up[n-1/2] <- up*[n-1/2] + dt/2 d/dz visc d/dz up[n-1/2]
  call cpu_clock_begin(id_clock_vertvisc)
  call enable_averages(dt, time_local, cs%diag)
  dt_visc = dt_pred ; if (cs%use_correct_dt_visc) dt_visc = dt
  call set_viscous_ml(u_in, v_in, h_av, tv, forces, visc, dt_visc, g, gv, us, cs%set_visc_CSp)
  call disable_averaging(cs%diag)

  call vertvisc_coef(up, vp, h_av, forces, visc, dt_pred, g, gv, us, cs%vertvisc_CSp, cs%OBC)
  call vertvisc(up, vp, h_av, forces, visc, dt_pred, cs%OBC, cs%ADp, cs%CDp, &
                g, gv, us, cs%vertvisc_CSp)
  call cpu_clock_end(id_clock_vertvisc)
  call pass_vector(up, vp, g%Domain, clock=id_clock_pass)

! uh = up[n-1/2] * h[n-1/2]
! h_av = h + dt div . uh
  call cpu_clock_begin(id_clock_continuity)
  call continuity(up, vp, h_in, hp, uh, vh, dt, g, gv, us, cs%continuity_CSp, obc=cs%OBC)
  call cpu_clock_end(id_clock_continuity)
  call pass_var(hp, g%Domain, clock=id_clock_pass)
  call pass_vector(uh, vh, g%Domain, clock=id_clock_pass)

! h_av <- (h + hp)/2   (centered at n-1/2)
  do k=1,nz ; do j=js-2,je+2 ; do i=is-2,ie+2
    h_av(i,j,k) = (h_in(i,j,k) + hp(i,j,k)) * 0.5
  enddo ; enddo ; enddo

  if (cs%debug) &
    call mom_state_chksum("Predictor 1", up, vp, h_av, uh, vh, g, gv, us)

! CAu = -(f+zeta(up))/h_av vh + d/dx KE(up)  (function of up[n-1/2], h[n-1/2])
  call cpu_clock_begin(id_clock_cor)
  call coradcalc(up, vp, h_av, uh, vh, cs%CAu, cs%CAv, cs%OBC, cs%ADp, &
                 g, gv, us, cs%CoriolisAdv_CSp)
  call cpu_clock_end(id_clock_cor)
  if (associated(cs%OBC)) then
    call open_boundary_zero_normal_flow(cs%OBC, g, cs%CAu, cs%CAv)
  endif

! call enable_averages(dt, Time_local, CS%diag)  ?????????????????????/

! up* = u[n] + (1+gamma) * dt * ( PFu + CAu )  Extrapolated for damping
! u*[n+1] = u[n] + dt * ( PFu + CAu )
  do k=1,nz ; do j=js,je ; do i=isq,ieq
    up(i,j,k) = g%mask2dCu(i,j) * (u_in(i,j,k) + dt * (1.+cs%begw) * &
            ((cs%PFu(i,j,k) + cs%CAu(i,j,k)) + cs%diffu(i,j,k)))
    u_in(i,j,k) = g%mask2dCu(i,j) * (u_in(i,j,k) +  dt * &
            ((cs%PFu(i,j,k) + cs%CAu(i,j,k)) + cs%diffu(i,j,k)))
  enddo ; enddo ; enddo
  do k=1,nz ; do j=jsq,jeq ; do i=is,ie
    vp(i,j,k) = g%mask2dCv(i,j) * (v_in(i,j,k) + dt * (1.+cs%begw) * &
            ((cs%PFv(i,j,k) + cs%CAv(i,j,k)) + cs%diffv(i,j,k)))
    v_in(i,j,k) = g%mask2dCv(i,j) * (v_in(i,j,k) + dt * &
            ((cs%PFv(i,j,k) + cs%CAv(i,j,k)) + cs%diffv(i,j,k)))
  enddo ; enddo ; enddo

! up[n] <- up* + dt d/dz visc d/dz up
! u[n] <- u*[n] + dt d/dz visc d/dz u[n]
  call cpu_clock_begin(id_clock_vertvisc)
  call vertvisc_coef(up, vp, h_av, forces, visc, dt, g, gv, us, cs%vertvisc_CSp, cs%OBC)
  call vertvisc(up, vp, h_av, forces, visc, dt, cs%OBC, cs%ADp, cs%CDp, &
                g, gv, us, cs%vertvisc_CSp, cs%taux_bot, cs%tauy_bot)
  call vertvisc_coef(u_in, v_in, h_av, forces, visc, dt, g, gv, us, cs%vertvisc_CSp, cs%OBC)
  call vertvisc(u_in, v_in, h_av, forces, visc, dt, cs%OBC, cs%ADp, cs%CDp,&
                g, gv, us, cs%vertvisc_CSp, cs%taux_bot, cs%tauy_bot)
  call cpu_clock_end(id_clock_vertvisc)
  call pass_vector(up, vp, g%Domain, clock=id_clock_pass)
  call pass_vector(u_in, v_in, g%Domain, clock=id_clock_pass)

! uh = up[n] * h[n]  (up[n] might be extrapolated to damp GWs)
! h[n+1] = h[n] + dt div . uh
  call cpu_clock_begin(id_clock_continuity)
  call continuity(up, vp, h_in, h_in, uh, vh, dt, g, gv, us, cs%continuity_CSp, obc=cs%OBC)
  call cpu_clock_end(id_clock_continuity)
  call pass_var(h_in, g%Domain, clock=id_clock_pass)
  call pass_vector(uh, vh, g%Domain, clock=id_clock_pass)

! Accumulate mass flux for tracer transport
  do k=1,nz
    do j=js-2,je+2 ; do i=isq-2,ieq+2
      uhtr(i,j,k) = uhtr(i,j,k) + dt*uh(i,j,k)
    enddo ; enddo
    do j=jsq-2,jeq+2 ; do i=is-2,ie+2
      vhtr(i,j,k) = vhtr(i,j,k) + dt*vh(i,j,k)
    enddo ; enddo
  enddo

  if (cs%debug) then
    call mom_state_chksum("Corrector", u_in, v_in, h_in, uh, vh, g, gv, us)
    call mom_accel_chksum("Corrector accel", cs%CAu, cs%CAv, cs%PFu, cs%PFv, &
                          cs%diffu, cs%diffv, g, gv, us)
  endif

  if (gv%Boussinesq) then
    do j=js,je ; do i=is,ie ; eta_av(i,j) = -gv%Z_to_H*g%bathyT(i,j) ; enddo ; enddo
  else
    do j=js,je ; do i=is,ie ; eta_av(i,j) = 0.0 ; enddo ; enddo
  endif
  do k=1,nz ; do j=js,je ; do i=is,ie
    eta_av(i,j) = eta_av(i,j) + h_av(i,j,k)
  enddo ; enddo ; enddo

  if (dyn_p_surf) deallocate(p_surf)

!   Here various terms used in to update the momentum equations are
! offered for averaging.
  if (cs%id_PFu > 0) call post_data(cs%id_PFu, cs%PFu, cs%diag)
  if (cs%id_PFv > 0) call post_data(cs%id_PFv, cs%PFv, cs%diag)
  if (cs%id_CAu > 0) call post_data(cs%id_CAu, cs%CAu, cs%diag)
  if (cs%id_CAv > 0) call post_data(cs%id_CAv, cs%CAv, cs%diag)

!   Here the thickness fluxes are offered for averaging.
  if (cs%id_uh > 0) call post_data(cs%id_uh, uh, cs%diag)
  if (cs%id_vh > 0) call post_data(cs%id_vh, vh, cs%diag)

end subroutine step_mom_dyn_unsplit_rk2

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

!> Allocate the control structure for this module, allocates memory in it, and registers
!! any auxiliary restart variables that are specific to the unsplit RK2 time stepping scheme.
!!
!! All variables registered here should have the ability to be recreated if they are not present
!! in a restart file.
subroutine register_restarts_dyn_unsplit_rk2(HI, GV, param_file, CS, restart_CS)
  type(hor_index_type),         intent(in)    :: HI         !< A horizontal index type structure.
  type(verticalgrid_type),      intent(in)    :: GV         !< The ocean's vertical grid structure.
  type(param_file_type),        intent(in)    :: param_file !< A structure to parse for run-time
                                                            !! parameters.
  type(mom_dyn_unsplit_rk2_cs), pointer       :: CS         !< The control structure set up by
                                                            !! initialize_dyn_unsplit_RK2.
  type(mom_restart_cs),         pointer       :: restart_CS !< A pointer to the restart control
                                                            !! structure.
!   This subroutine sets up any auxiliary restart variables that are specific
! to the unsplit time stepping scheme.  All variables registered here should
! have the ability to be recreated if they are not present in a restart file.

  ! Local variables
  character(len=48) :: thickness_units, flux_units
  integer :: isd, ied, jsd, jed, nz, IsdB, IedB, JsdB, JedB
  isd = hi%isd ; ied = hi%ied ; jsd = hi%jsd ; jed = hi%jed ; nz = gv%ke
  isdb = hi%IsdB ; iedb = hi%IedB ; jsdb = hi%JsdB ; jedb = hi%JedB

! This is where a control structure that is specific to this module would be allocated.
  if (associated(cs)) then
    call mom_error(warning, "register_restarts_dyn_unsplit_RK2 called with an associated "// &
                             "control structure.")
    return
  endif
  allocate(cs)

  alloc_(cs%diffu(isdb:iedb,jsd:jed,nz)) ; cs%diffu(:,:,:) = 0.0
  alloc_(cs%diffv(isd:ied,jsdb:jedb,nz)) ; cs%diffv(:,:,:) = 0.0
  alloc_(cs%CAu(isdb:iedb,jsd:jed,nz)) ; cs%CAu(:,:,:) = 0.0
  alloc_(cs%CAv(isd:ied,jsdb:jedb,nz)) ; cs%CAv(:,:,:) = 0.0
  alloc_(cs%PFu(isdb:iedb,jsd:jed,nz)) ; cs%PFu(:,:,:) = 0.0
  alloc_(cs%PFv(isd:ied,jsdb:jedb,nz)) ; cs%PFv(:,:,:) = 0.0

  thickness_units = get_thickness_units(gv)
  flux_units = get_flux_units(gv)

!  No extra restart fields are needed with this time stepping scheme.

end subroutine register_restarts_dyn_unsplit_rk2

!> Initialize parameters and allocate memory associated with the unsplit RK2 dynamics module.
subroutine initialize_dyn_unsplit_rk2(u, v, h, Time, G, GV, US, param_file, diag, CS, &
                                      restart_CS, Accel_diag, Cont_diag, MIS, MEKE, &
                                      OBC, update_OBC_CSp, ALE_CSp, setVisc_CSp, &
                                      visc, dirs, ntrunc, cont_stencil)
  type(ocean_grid_type),                     intent(inout) :: G    !< The ocean's grid structure.
  type(verticalgrid_type),                   intent(in)    :: GV   !< The ocean's vertical grid structure.
  type(unit_scale_type),                     intent(in)    :: US   !< A dimensional unit scaling type
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: u    !< The zonal velocity [L T-1 ~> m s-1].
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: v    !< The meridional velocity [L T-1 ~> m s-1].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) , intent(inout) :: h    !< Layer thicknesses [H ~> m or kg m-2]
  type(time_type),                   target, intent(in)    :: Time !< The current model time.
  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(mom_dyn_unsplit_rk2_cs),              pointer       :: CS   !< The control structure set up
                                                                   !! by initialize_dyn_unsplit_RK2.
  type(mom_restart_cs),                      pointer       :: restart_CS !< A pointer to the restart
                                                                         !! control structure.
  type(accel_diag_ptrs),             target, intent(inout) :: Accel_diag !< A set of pointers to the
                                      !! various accelerations in the momentum equations, which can
                                      !! be used for later derived diagnostics, like energy budgets.
  type(cont_diag_ptrs),              target, intent(inout) :: Cont_diag  !< A structure with pointers
                                                                         !! to various terms in the
                                                                         !! continuity equations.
  type(ocean_internal_state),                intent(inout) :: MIS  !< The "MOM6 Internal State"
                                                     !! structure, used to pass around pointers
                                                     !! to various arrays for diagnostic purposes.
  type(meke_type),                           pointer       :: MEKE !< MEKE data
  type(ocean_obc_type),                      pointer       :: OBC  !< If open boundary conditions
                                                    !! are used, this points to the ocean_OBC_type
                                                    !! that was set up in MOM_initialization.
  type(update_obc_cs),                       pointer       :: update_OBC_CSp !< If open boundary
                                                         !! condition updates are used, this points
                                                         !! to the appropriate control structure.
  type(ale_cs),                              pointer       :: ALE_CSp     !< This points to the ALE
                                                                          !! control structure.
  type(set_visc_cs),                         pointer       :: setVisc_CSp !< This points to the
                                                                          !! set_visc control
                                                                          !! structure.
  type(vertvisc_type),                       intent(inout) :: visc !< A structure containing
                                                         !! vertical viscosities, bottom drag
                                                         !! viscosities, and related fields.
  type(directories),                         intent(in)    :: dirs !< A structure containing several
                                                                   !! relevant directory paths.
  integer, target,                           intent(inout) :: ntrunc !< A target for the variable
                                                       !! that records the number of times the
                                                       !! velocity is truncated (this should be 0).
  integer,                         optional, intent(out)   :: cont_stencil !< The stencil for
                                                       !! thickness from the continuity solver.

  !   This subroutine initializes all of the variables that are used by this
  ! dynamic core, including diagnostics and the cpu clocks.

  ! Local variables
  character(len=40) :: mdl = "MOM_dynamics_unsplit_RK2" ! This module's name.
  character(len=48) :: thickness_units, flux_units
  ! This include declares and sets the variable "version".
# include "version_variable.h"
  real :: H_convert
  logical :: use_tides
  integer :: isd, ied, jsd, jed, nz, IsdB, IedB, JsdB, JedB
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = g%ke
  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB

  if (.not.associated(cs)) call mom_error(fatal, &
      "initialize_dyn_unsplit_RK2 called with an unassociated control structure.")
  if (cs%module_is_initialized) then
    call mom_error(warning, "initialize_dyn_unsplit_RK2 called with a control "// &
                            "structure that has already been initialized.")
    return
  endif
  cs%module_is_initialized = .true.

  cs%diag => diag

  call log_version(param_file, mdl, version, "")
  call get_param(param_file, mdl, "BE", cs%be, &
                 "If SPLIT is true, BE determines the relative weighting "//&
                 "of a  2nd-order Runga-Kutta baroclinic time stepping "//&
                 "scheme (0.5) and a backward Euler scheme (1) that is "//&
                 "used for the Coriolis and inertial terms.  BE may be "//&
                 "from 0.5 to 1, but instability may occur near 0.5. "//&
                 "BE is also applicable if SPLIT is false and USE_RK2 "//&
                 "is true.", units="nondim", default=0.6)
  call get_param(param_file, mdl, "BEGW", cs%begw, &
                 "If SPLIT is true, BEGW is a number from 0 to 1 that "//&
                 "controls the extent to which the treatment of gravity "//&
                 "waves is forward-backward (0) or simulated backward "//&
                 "Euler (1).  0 is almost always used. "//&
                 "If SPLIT is false and USE_RK2 is true, BEGW can be "//&
                 "between 0 and 0.5 to damp gravity waves.", &
                 units="nondim", default=0.0)
  call get_param(param_file, mdl, "FIX_UNSPLIT_DT_VISC_BUG", cs%use_correct_dt_visc, &
                 "If true, use the correct timestep in the viscous terms applied in the first "//&
                 "predictor step with the unsplit time stepping scheme, and in the calculation "//&
                 "of the turbulent mixed layer properties for viscosity with unsplit or "//&
                 "unsplit_RK2.", default=.true.)
  call get_param(param_file, mdl, "DEBUG", cs%debug, &
                 "If true, write out verbose debugging data.", &
                 default=.false., debuggingparam=.true.)
  call get_param(param_file, mdl, "TIDES", use_tides, &
                 "If true, apply tidal momentum forcing.", default=.false.)

  allocate(cs%taux_bot(isdb:iedb,jsd:jed)) ; cs%taux_bot(:,:) = 0.0
  allocate(cs%tauy_bot(isd:ied,jsdb:jedb)) ; cs%tauy_bot(:,:) = 0.0

  mis%diffu => cs%diffu ; mis%diffv => cs%diffv
  mis%PFu => cs%PFu ; mis%PFv => cs%PFv
  mis%CAu => cs%CAu ; mis%CAv => cs%CAv

  cs%ADp => accel_diag ; cs%CDp => cont_diag
  accel_diag%diffu => cs%diffu ; accel_diag%diffv => cs%diffv
  accel_diag%PFu => cs%PFu ; accel_diag%PFv => cs%PFv
  accel_diag%CAu => cs%CAu ; accel_diag%CAv => cs%CAv

  call continuity_init(time, g, gv, us, param_file, diag, cs%continuity_CSp)
  if (present(cont_stencil)) cont_stencil = continuity_stencil(cs%continuity_CSp)
  call coriolisadv_init(time, g, gv, us, param_file, diag, cs%ADp, cs%CoriolisAdv_CSp)
  if (use_tides) call tidal_forcing_init(time, g, param_file, cs%tides_CSp)
  call pressureforce_init(time, g, gv, us, param_file, diag, cs%PressureForce_CSp, &
                          cs%tides_CSp)
  call hor_visc_init(time, g, us, param_file, diag, cs%hor_visc_CSp, meke)
  call vertvisc_init(mis, time, g, gv, us, param_file, diag, cs%ADp, dirs, &
                     ntrunc, cs%vertvisc_CSp)
  if (.not.associated(setvisc_csp)) call mom_error(fatal, &
    "initialize_dyn_unsplit_RK2 called with setVisc_CSp unassociated.")
  cs%set_visc_CSp => setvisc_csp

  if (associated(ale_csp)) cs%ALE_CSp => ale_csp
  if (associated(obc)) cs%OBC => obc

  flux_units = get_flux_units(gv)
  h_convert = gv%H_to_m ; if (.not.gv%Boussinesq) h_convert = gv%H_to_kg_m2
  cs%id_uh = register_diag_field('ocean_model', 'uh', diag%axesCuL, time, &
      'Zonal Thickness Flux', flux_units, y_cell_method='sum', v_extensive=.true., &
      conversion=h_convert*us%L_to_m**2*us%s_to_T)
  cs%id_vh = register_diag_field('ocean_model', 'vh', diag%axesCvL, time, &
      'Meridional Thickness Flux', flux_units, x_cell_method='sum', v_extensive=.true., &
      conversion=h_convert*us%L_to_m**2*us%s_to_T)
  cs%id_CAu = register_diag_field('ocean_model', 'CAu', diag%axesCuL, time, &
      'Zonal Coriolis and Advective Acceleration', 'meter second-2', conversion=us%L_T2_to_m_s2)
  cs%id_CAv = register_diag_field('ocean_model', 'CAv', diag%axesCvL, time, &
      'Meridional Coriolis and Advective Acceleration', 'meter second-2', conversion=us%L_T2_to_m_s2)
  cs%id_PFu = register_diag_field('ocean_model', 'PFu', diag%axesCuL, time, &
      'Zonal Pressure Force Acceleration', 'meter second-2', conversion=us%L_T2_to_m_s2)
  cs%id_PFv = register_diag_field('ocean_model', 'PFv', diag%axesCvL, time, &
      'Meridional Pressure Force Acceleration', 'meter second-2', conversion=us%L_T2_to_m_s2)

  id_clock_cor = cpu_clock_id('(Ocean Coriolis & mom advection)', grain=clock_module)
  id_clock_continuity = cpu_clock_id('(Ocean continuity equation)', grain=clock_module)
  id_clock_pres = cpu_clock_id('(Ocean pressure force)', grain=clock_module)
  id_clock_vertvisc = cpu_clock_id('(Ocean vertical viscosity)', grain=clock_module)
  id_clock_horvisc = cpu_clock_id('(Ocean horizontal viscosity)', grain=clock_module)
  id_clock_mom_update = cpu_clock_id('(Ocean momentum increments)', grain=clock_module)
  id_clock_pass = cpu_clock_id('(Ocean message passing)', grain=clock_module)
  id_clock_pass_init = cpu_clock_id('(Ocean init message passing)', grain=clock_routine)

end subroutine initialize_dyn_unsplit_rk2

!> Clean up and deallocate memory associated with the dyn_unsplit_RK2 module.
subroutine end_dyn_unsplit_rk2(CS)
  type(mom_dyn_unsplit_rk2_cs), pointer :: CS !< dyn_unsplit_RK2 control structure that
                                              !! will be deallocated in this subroutine.

  dealloc_(cs%diffu) ; dealloc_(cs%diffv)
  dealloc_(cs%CAu)   ; dealloc_(cs%CAv)
  dealloc_(cs%PFu)   ; dealloc_(cs%PFv)

  deallocate(cs)
end subroutine end_dyn_unsplit_rk2

end module mom_dynamics_unsplit_rk2