MOM_dynamics_unsplit.F90

!> Time steps the ocean dynamics with an unsplit quasi 3rd order scheme
module mom_dynamics_unsplit

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

!********+*********+*********+*********+*********+*********+*********+**
!*                                                                     *
!*  By Robert Hallberg, 1993-2012                                      *
!*                                                                     *
!*    This file contains code that does the time-stepping of the       *
!*  adiabatic dynamic core, in this case with an unsplit third-order   *
!*  Runge-Kutta time stepping scheme for the momentum and a forward-   *
!*  backward coupling between the momentum and continuity equations.   *
!*  This was the orignal unsplit time stepping scheme used in early    *
!*  versions of HIM and its precuror.  While it is very simple and     *
!*  accurate, it is much less efficient that the split time stepping   *
!*  scheme for realistic oceanographic applications.  It has been      *
!*  retained for all of these years primarily to verify that the split *
!*  scheme is giving the right answers, and to debug the failings of   *
!*  the split scheme when it is not.  The split time stepping scheme   *
!*  is now sufficiently robust that it should be first choice for      *
!*  almost any conceivable application, except perhaps from cases      *
!*  with just a few layers for which the exact timing of the high-     *
!*  frequency barotropic gravity waves is of paramount importance.     *
!*  This scheme is slightly more efficient than the other unsplit      *
!*  scheme that can be found in MOM_dynamics_unsplit_RK2.F90.          *
!*                                                                     *
!*    The subroutine step_MOM_dyn_unsplit actually does the time       *
!*  stepping, while register_restarts_dyn_unsplit  sets the fields     *
!*  that are found in a full restart file with this scheme, and        *
!*  initialize_dyn_unsplit  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 : accel_diag_ptrs, ocean_internal_state, 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, diag_update_remap_grids
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_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, real_to_time, operator(+)
use mom_time_manager, only : operator(-), operator(>), operator(*), operator(/)

use mom_ale, only : ale_cs
use mom_barotropic, only : barotropic_cs
use mom_boundary_update, only : update_obc_data, update_obc_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_interface_heights, only : find_eta
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
use mom_wave_interface, only: wave_parameters_cs

implicit none ; private

#include <MOM_memory.h>

!> MOM_dynamics_unsplit module control structure
type, public :: mom_dyn_unsplit_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]

  logical :: 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.  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_cs

public step_mom_dyn_unsplit, register_restarts_dyn_unsplit
public initialize_dyn_unsplit, end_dyn_unsplit

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

contains

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

!> Step the MOM6 dynamics using an unsplit mixed 2nd order (for continuity) and
!! 3rd order (for the inviscid momentum equations) order scheme
subroutine step_mom_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, &
                  p_surf_begin, p_surf_end, uh, vh, uhtr, vhtr, eta_av, G, GV, US, CS, &
                  VarMix, MEKE, Waves)
  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(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 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 start 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_cs), pointer      :: CS     !< The control structure set up by
                                                   !! initialize_dyn_unsplit.
  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.
  type(wave_parameters_cs), optional, pointer :: Waves !< A pointer to a structure containing
                                 !! fields related to the surface wave conditions

  ! Local variables
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_av, hp ! Prediced or averaged layer thicknesses [H ~> m or kg m-2]
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: up, upp ! Predicted zonal velocities [L T-1 ~> m s-1]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: vp, vpp ! 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 / 3.0

  h_av(:,:,:) = 0; hp(:,:,:) = 0
  up(:,:,:) = 0; upp(:,:,:) = 0
  vp(:,:,:) = 0; vpp(:,:,:) = 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

! Matsuno's third order accurate three 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 First Predictor ", u, v, h, 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, v, h, 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)

! uh = u*h
! hp = h + dt/2 div . uh
  call cpu_clock_begin(id_clock_continuity)
  call continuity(u, v, h, hp, uh, vh, dt*0.5, 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)

  call enable_averages(0.5*dt, time_local-real_to_time(0.5*us%T_to_s*dt), cs%diag)
!   Here the first half of 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)
  call disable_averaging(cs%diag)

! h_av = (h + hp)/2
! u = u + dt diffu
  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(i,j,k) + hp(i,j,k)) * 0.5
    enddo ; enddo
    do j=js,je ; do i=isq,ieq
      u(i,j,k) = u(i,j,k) + dt * cs%diffu(i,j,k) * g%mask2dCu(i,j)
    enddo ; enddo
    do j=jsq,jeq ; do i=is,ie
      v(i,j,k) = v(i,j,k) + dt * cs%diffv(i,j,k) * g%mask2dCv(i,j)
    enddo ; enddo
    do j=js-2,je+2 ; do i=isq-2,ieq+2
      uhtr(i,j,k) = uhtr(i,j,k) + 0.5*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) + 0.5*dt*vh(i,j,k)
    enddo ; enddo
  enddo
  call cpu_clock_end(id_clock_mom_update)
  call pass_vector(u, v, g%Domain, clock=id_clock_pass)

! CAu = -(f+zeta)/h_av vh + d/dx KE
  call cpu_clock_begin(id_clock_cor)
  call coradcalc(u, v, 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)
  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.75*p_surf_begin(i,j) + 0.25*p_surf_end(i,j)
  enddo ; enddo ; endif
  call pressureforce(h_av, tv, cs%PFu, cs%PFv, g, gv, us, &
                     cs%PressureForce_CSp, cs%ALE_CSp, p_surf)
  call cpu_clock_end(id_clock_pres)

  if (associated(cs%OBC)) then; if (cs%OBC%update_OBC) then
    call update_obc_data(cs%OBC, g, gv, us, tv, h, 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)
  endif

! up = u + 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(i,j,k) + dt_pred * (cs%PFu(i,j,k) + cs%CAu(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(i,j,k) + dt_pred * (cs%PFv(i,j,k) + cs%CAv(i,j,k)))
  enddo ; enddo ; enddo
  call cpu_clock_end(id_clock_mom_update)

  if (cs%debug) then
    call mom_state_chksum("Predictor 1", up, vp, h_av, uh, vh, g, gv, us)
    call mom_accel_chksum("Predictor 1 accel", cs%CAu, cs%CAv, cs%PFu, cs%PFv, &
                          cs%diffu, cs%diffv, g, gv, us)
  endif

 ! up <- up + dt/2 d/dz visc d/dz up
  call cpu_clock_begin(id_clock_vertvisc)
  call enable_averages(dt, time_local, cs%diag)
  dt_visc = 0.5*dt ; if (cs%use_correct_dt_visc) dt_visc = dt
  call set_viscous_ml(u, v, h_av, tv, forces, visc, dt_visc, g, gv, us, cs%set_visc_CSp)
  call disable_averaging(cs%diag)

  dt_visc = 0.5*dt ; if (cs%use_correct_dt_visc) dt_visc = dt_pred
  call vertvisc_coef(up, vp, h_av, forces, visc, dt_visc, g, gv, us, cs%vertvisc_CSp, cs%OBC)
  call vertvisc(up, vp, h_av, forces, visc, dt_visc, cs%OBC, cs%ADp, cs%CDp, &
                g, gv, us, cs%vertvisc_CSp, waves=waves)
  call cpu_clock_end(id_clock_vertvisc)
  call pass_vector(up, vp, g%Domain, clock=id_clock_pass)

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

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

! CAu = -(f+zeta(up))/h_av vh + d/dx KE(up)
  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)

! PFu = d/dx M(h_av,T,S)
  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.25*p_surf_begin(i,j) + 0.75*p_surf_end(i,j)
  enddo ; enddo ; endif
  call pressureforce(h_av, tv, cs%PFu, cs%PFv, g, gv, us, &
                     cs%PressureForce_CSp, cs%ALE_CSp, p_surf)
  call cpu_clock_end(id_clock_pres)

  if (associated(cs%OBC)) then; if (cs%OBC%update_OBC) then
    call update_obc_data(cs%OBC, g, gv, us, tv, h, 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)
  endif

! upp = u + dt/2 * ( PFu + CAu )
  call cpu_clock_begin(id_clock_mom_update)
  do k=1,nz ; do j=js,je ; do i=isq,ieq
    upp(i,j,k) = g%mask2dCu(i,j) * (u(i,j,k) + dt * 0.5 * (cs%PFu(i,j,k) + cs%CAu(i,j,k)))
  enddo ; enddo ; enddo
  do k=1,nz ; do j=jsq,jeq ; do i=is,ie
    vpp(i,j,k) = g%mask2dCv(i,j) * (v(i,j,k) + dt * 0.5 * (cs%PFv(i,j,k) + cs%CAv(i,j,k)))
  enddo ; enddo ; enddo
  call cpu_clock_end(id_clock_mom_update)

  if (cs%debug) then
    call mom_state_chksum("Predictor 2", upp, vpp, h_av, uh, vh, g, gv, us)
    call mom_accel_chksum("Predictor 2 accel", cs%CAu, cs%CAv, cs%PFu, cs%PFv, &
                          cs%diffu, cs%diffv, g, gv, us)
  endif

! upp <- upp + dt/2 d/dz visc d/dz upp
  call cpu_clock_begin(id_clock_vertvisc)
  call vertvisc_coef(upp, vpp, hp, forces, visc, dt*0.5, g, gv, us, cs%vertvisc_CSp, cs%OBC)
  call vertvisc(upp, vpp, hp, forces, visc, dt*0.5, cs%OBC, cs%ADp, cs%CDp, &
                g, gv, us, cs%vertvisc_CSp, waves=waves)
  call cpu_clock_end(id_clock_vertvisc)
  call pass_vector(upp, vpp, g%Domain, clock=id_clock_pass)

! uh = upp * hp
! h = hp + dt/2 div . uh
  call cpu_clock_begin(id_clock_continuity)
  call continuity(upp, vpp, hp, h, uh, vh, (dt*0.5), g, gv, us, cs%continuity_CSp, obc=cs%OBC)
  call cpu_clock_end(id_clock_continuity)
  call pass_var(h, g%Domain, clock=id_clock_pass)
  call pass_vector(uh, vh, g%Domain, clock=id_clock_pass)
  ! Whenever thickness changes let the diag manager know, target grids
  ! for vertical remapping may need to be regenerated.
  call diag_update_remap_grids(cs%diag)

  call enable_averages(0.5*dt, time_local, cs%diag)
!   Here the second half of 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)
  call disable_averaging(cs%diag)
  call enable_averages(dt, time_local, cs%diag)

! h_av = (h + hp)/2
  do k=1,nz
    do j=js-2,je+2 ; do i=is-2,ie+2
      h_av(i,j,k) = 0.5*(h(i,j,k) + hp(i,j,k))
    enddo ; enddo
    do j=js-2,je+2 ; do i=isq-2,ieq+2
      uhtr(i,j,k) = uhtr(i,j,k) + 0.5*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) + 0.5*dt*vh(i,j,k)
    enddo ; enddo
  enddo

! CAu = -(f+zeta(upp))/h_av vh + d/dx KE(upp)
  call cpu_clock_begin(id_clock_cor)
  call coradcalc(upp, vpp, 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)
  call cpu_clock_begin(id_clock_pres)
  call pressureforce(h_av, tv, cs%PFu, cs%PFv, g, gv, us, &
                     cs%PressureForce_CSp, cs%ALE_CSp, p_surf)
  call cpu_clock_end(id_clock_pres)

  if (associated(cs%OBC)) then; if (cs%OBC%update_OBC) then
    call update_obc_data(cs%OBC, g, gv, us, tv, h, cs%update_OBC_CSp, time_local)
  endif; endif

! u = u + dt * ( PFu + CAu )
  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)
  endif
  do k=1,nz ; do j=js,je ; do i=isq,ieq
    u(i,j,k) = g%mask2dCu(i,j) * (u(i,j,k) + dt * (cs%PFu(i,j,k) + cs%CAu(i,j,k)))
  enddo ; enddo ; enddo
  do k=1,nz ; do j=jsq,jeq ; do i=is,ie
    v(i,j,k) = g%mask2dCv(i,j) * (v(i,j,k) + dt * (cs%PFv(i,j,k) + cs%CAv(i,j,k)))
  enddo ; enddo ; enddo

! u <- u + dt d/dz visc d/dz u
  call cpu_clock_begin(id_clock_vertvisc)
  call vertvisc_coef(u, v, h_av, forces, visc, dt, g, gv, us, cs%vertvisc_CSp, cs%OBC)
  call vertvisc(u, v, h_av, forces, visc, dt, cs%OBC, cs%ADp, cs%CDp, &
                g, gv, us, cs%vertvisc_CSp, cs%taux_bot, cs%tauy_bot, waves=waves)
  call cpu_clock_end(id_clock_vertvisc)
  call pass_vector(u, v, g%Domain, clock=id_clock_pass)

  if (cs%debug) then
    call mom_state_chksum("Corrector", u, v, h, 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)

end subroutine step_mom_dyn_unsplit

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

!> Allocate the control structure for this module, allocates memory in it, and registers
!! 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.
subroutine register_restarts_dyn_unsplit(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_cs),  pointer    :: CS         !< The control structure set up by
                                                      !! initialize_dyn_unsplit.
  type(mom_restart_cs),      pointer    :: restart_CS !< A pointer to the restart control structure.

  ! Local arguments
  character(len=40)  :: mdl = "MOM_dynamics_unsplit" ! This module's name.
  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 is allocated.
  if (associated(cs)) then
    call mom_error(warning, "register_restarts_dyn_unsplit 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

!> Initialize parameters and allocate memory associated with the unsplit dynamics module.
subroutine initialize_dyn_unsplit(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_cs),       pointer       :: CS         !< The control structure set up
                                                              !! by initialize_dyn_unsplit.
  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" ! 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 called with an unassociated control structure.")
  if (cs%module_is_initialized) then
    call mom_error(warning, "initialize_dyn_unsplit 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, "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 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
  if (associated(update_obc_csp)) cs%update_OBC_CSp => update_obc_csp

  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', 'm s-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', 'm s-2', &
      conversion=us%L_T2_to_m_s2)
  cs%id_PFu = register_diag_field('ocean_model', 'PFu', diag%axesCuL, time, &
      'Zonal Pressure Force Acceleration', 'm s-2', &
      conversion=us%L_T2_to_m_s2)
  cs%id_PFv = register_diag_field('ocean_model', 'PFv', diag%axesCvL, time, &
      'Meridional Pressure Force Acceleration', 'm s-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

!> Clean up and deallocate memory associated with the unsplit dynamics module.
subroutine end_dyn_unsplit(CS)
  type(mom_dyn_unsplit_cs), pointer :: CS !< unsplit dynamics 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

end module mom_dynamics_unsplit