MOM_dynamics_unsplit.F90

Go to the documentation of this file.

module mom_dynamics_unsplit

!***********************************************************************
!*                   GNU General Public License                        *
!* This file is a part of MOM.                                         *
!*                                                                     *
!* MOM is free software; you can redistribute it and/or modify it and  *
!* are expected to follow the terms of the GNU General Public License  *
!* as published by the Free Software Foundation; either version 2 of   *
!* the License, or (at your option) any later version.                 *
!*                                                                     *
!* MOM is distributed in the hope that it will be useful, but WITHOUT  *
!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY  *
!* or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public    *
!* License for more details.                                           *
!*                                                                     *
!* For the full text of the GNU General Public License,                *
!* write to: Free Software Foundation, Inc.,                           *
!*           675 Mass Ave, Cambridge, MA 02139, USA.                   *
!* or see:   http://www.gnu.org/licenses/gpl.html                      *
!***********************************************************************

!********+*********+*********+*********+*********+*********+*********+**
!*                                                                     *
!*  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 : 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_averaging
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, vardesc
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, set_time, 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_continuity, only : continuity, continuity_init, continuity_cs
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_tidal_forcing, only : tidal_forcing_init, tidal_forcing_cs
use mom_vert_friction, only : vertvisc, vertvisc_coef
use mom_vert_friction, only : vertvisc_limit_vel, 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>
type, public :: mom_dyn_unsplit_cs ; private
  real allocable_, dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: &
    cau, &    ! CAu = f*v - u.grad(u) in m s-2.
    pfu, &    ! PFu = -dM/dx, in m s-2.
    diffu     ! Zonal acceleration due to convergence of the along-isopycnal
              ! stress tensor, in m s-2.

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

  real, pointer, dimension(:,:) :: taux_bot => null(), tauy_bot => null()
    ! The frictional bottom stresses from the ocean to the seafloor, in Pa.

  logical :: debug           ! If true, write verbose checksums for debugging purposes.

  logical :: module_is_initialized = .false.

  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 ! A structure that is used to regulate the
                                   ! timing of diagnostic output.
  type(accel_diag_ptrs), pointer :: adp ! A structure pointing to the various
                                   ! accelerations in the momentum equations,
                                   ! which can later be used to calculate
                                   ! derived diagnostics like energy budgets.
  type(cont_diag_ptrs), pointer :: cdp ! 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 is pointers to child subroutines' control strings.
  type(hor_visc_cs), pointer :: hor_visc_csp => null()
  type(continuity_cs), pointer :: continuity_csp => null()
  type(coriolisadv_cs), pointer :: coriolisadv_csp => null()
  type(pressureforce_cs), pointer :: pressureforce_csp => null()
  type(vertvisc_cs), pointer :: vertvisc_csp => null()
  type(set_visc_cs), pointer :: set_visc_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.
  type(update_obc_cs),    pointer :: update_obc_csp => null()
  type(tidal_forcing_cs), pointer :: tides_csp => null()

! This is a copy of the pointer in the top-level control structure.
  type(ale_cs), pointer :: ale_csp => null()

end type mom_dyn_unsplit_cs

public step_mom_dyn_unsplit, register_restarts_dyn_unsplit
public initialize_dyn_unsplit, end_dyn_unsplit

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

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

subroutine step_mom_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, fluxes, &
                  p_surf_begin, p_surf_end, uh, vh, uhtr, vhtr, eta_av, G, GV, 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.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                                    intent(inout) :: u      !< The zonal velocity, in m s-1.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                                    intent(inout) :: v      !< The meridional velocity, in m s-1.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  &
                                    intent(inout) :: h      !< Layer thicknesses, in H.
                                                            !! (usually 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, in s.
  type(forcing),                    intent(in)    :: fluxes !< A structure containing pointers to
                                 !! any possible forcing fields.  Unused fields have NULL ptrs.
  real, dimension(:,:),             pointer       :: p_surf_begin !< A pointer (perhaps NULL) to the
                                 !! surface pressure at the beginning of this dynamic step, in Pa.
  real, dimension(:,:),             pointer       :: p_surf_end   !< A pointer (perhaps NULL) to the
                                 !! surface pressure at the end of this dynamic step, in Pa.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                                    intent(inout) :: uh     !< The zonal volume or mass transport,
                                                            !! in m3 s-1 or kg s-1.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                                    intent(inout) :: vh     !< The meridional volume or mass
                                                            !! transport, in m3 s-1 or kg s-1.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                                    intent(inout) :: uhtr   !< he accumulated zonal volume or mass
                                        !! transport since the last tracer advection, in 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, in m3 or kg.
  real, dimension(SZI_(G),SZJ_(G)), intent(out)   :: eta_av !< The time-mean free surface height or
                                                            !! column mass, in 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.
! Arguments: u - The input and output zonal velocity, in m s-1.
!  (inout)   v - The input and output meridional velocity, in m s-1.
!  (inout)   h - The input and output layer thicknesses, in m or kg m-2,
!                depending on whether the Boussinesq approximation is made.
!  (in)      tv - a structure pointing to various thermodynamic variables.
!  (inout)   visc - A structure containing vertical viscosities, bottom drag
!                   viscosities, and related fields.
!  (in)      Time_local - The model time at the end of the time step.
!  (in)      dt - The time step in s.
!  (in)      fluxes - A structure containing pointers to any possible
!                     forcing fields.  Unused fields have NULL ptrs.
!  (in)      p_surf_begin - A pointer (perhaps NULL) to the surface pressure
!                     at the beginning of this dynamic step, in Pa.
!  (in)      p_surf_end - A pointer (perhaps NULL) to the surface pressure
!                     at the end of this dynamic step, in Pa.
!  (inout)   uh - The zonal volume or mass transport, in m3 s-1 or kg s-1.
!  (inout)   vh - The meridional volume or mass transport, in m3 s-1 or kg s-1.
!  (inout)   uhtr - The accumulated zonal volume or mass transport since the last
!                   tracer advection, in m3 or kg.
!  (inout)   vhtr - The accumulated meridional volume or mass transport since the last
!                   tracer advection, in m3 or kg.
!  (out)     eta_av - The time-mean free surface height or column mass, in m or
!                     kg m-2.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      CS - The control structure set up by initialize_dyn_unsplit.
!  (in)      VarMix - A pointer to a structure with fields that specify the
!                     spatially variable viscosities.
!  (inout)   MEKE - A pointer to a structure containing fields related to
!                   the Mesoscale Eddy Kinetic Energy.

  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_av, hp
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: up, upp
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: vp, vpp
  real, dimension(:,:), pointer :: p_surf
  real :: dt_pred   ! The time step for the predictor part of the baroclinic
                    ! time stepping.
  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 => fluxes%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)
  endif

! diffu = horizontal viscosity terms (u,h)
  call enable_averaging(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, 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, cs%continuity_CSp, &
                  obc=cs%OBC)
  call cpu_clock_end(id_clock_continuity)
  call cpu_clock_begin(id_clock_pass)
  call pass_var(hp, g%Domain)
  call pass_vector(uh, vh, g%Domain)
  call cpu_clock_end(id_clock_pass)

  call enable_averaging(0.5*dt,time_local-set_time(int(0.5*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 cpu_clock_begin(id_clock_pass)
  call pass_vector(u, v, g%Domain)
  call cpu_clock_end(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, 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, &
                     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, 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)
    call mom_accel_chksum("Predictor 1 accel", cs%CAu, cs%CAv, cs%PFu, cs%PFv,&
                          cs%diffu, cs%diffv, g, gv)
  endif

 ! up <- up + dt/2 d/dz visc d/dz up
  call cpu_clock_begin(id_clock_vertvisc)
  call enable_averaging(dt, time_local, cs%diag)
  call set_viscous_ml(u, v, h_av, tv, fluxes, visc, dt*0.5, g, gv, &
                      cs%set_visc_CSp)
  call disable_averaging(cs%diag)
  call vertvisc_coef(up, vp, h_av, fluxes, visc, dt*0.5, g, gv, &
                     cs%vertvisc_CSp, cs%OBC)
  call vertvisc(up, vp, h_av, fluxes, visc, dt*0.5, cs%OBC, cs%ADp, cs%CDp, &
                g, gv, cs%vertvisc_CSp)
  call cpu_clock_end(id_clock_vertvisc)
  call cpu_clock_begin(id_clock_pass)
  call pass_vector(up, vp, g%Domain)
  call cpu_clock_end(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, cs%continuity_CSp, obc=cs%OBC)
  call cpu_clock_end(id_clock_continuity)
  call cpu_clock_begin(id_clock_pass)
  call pass_var(h_av, g%Domain)
  call pass_vector(uh, vh, g%Domain)
  call cpu_clock_end(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, 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, &
                     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, 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)
    call mom_accel_chksum("Predictor 2 accel", cs%CAu, cs%CAv, cs%PFu, cs%PFv,&
                          cs%diffu, cs%diffv, g, gv)
  endif

! upp <- upp + dt/2 d/dz visc d/dz upp
  call cpu_clock_begin(id_clock_vertvisc)
  call vertvisc_coef(upp, vpp, hp, fluxes, visc, dt*0.5, g, gv, &
                     cs%vertvisc_CSp, cs%OBC)
  call vertvisc(upp, vpp, hp, fluxes, visc, dt*0.5, cs%OBC, cs%ADp, cs%CDp, &
                g, gv, cs%vertvisc_CSp)
  call cpu_clock_end(id_clock_vertvisc)
  call cpu_clock_begin(id_clock_pass)
  call pass_vector(upp, vpp, g%Domain)
  call cpu_clock_end(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, cs%continuity_CSp, obc=cs%OBC)
  call cpu_clock_end(id_clock_continuity)
  call cpu_clock_begin(id_clock_pass)
  call pass_var(h, g%Domain)
  call pass_vector(uh, vh, g%Domain)
  call cpu_clock_end(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_averaging(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_averaging(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, 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, &
                     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, 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, fluxes, visc, dt, g, gv, cs%vertvisc_CSp, cs%OBC)
  call vertvisc(u, v, h_av, fluxes, visc, dt, cs%OBC, cs%ADp, cs%CDp, &
                g, gv, cs%vertvisc_CSp, cs%taux_bot, cs%tauy_bot)
  call cpu_clock_end(id_clock_vertvisc)
  call cpu_clock_begin(id_clock_pass)
  call pass_vector(u, v, g%Domain)
  call cpu_clock_end(id_clock_pass)

  if (cs%debug) then
    call mom_state_chksum("Corrector", u, v, h, uh, vh, g, gv)
    call mom_accel_chksum("Corrector accel", cs%CAu, cs%CAv, cs%PFu, cs%PFv, &
                          cs%diffu, cs%diffv, g, gv)
  endif

  if (gv%Boussinesq) then
    do j=js,je ; do i=is,ie ; eta_av(i,j) = -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

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

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.

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

! Arguments: HI - A horizontal index type structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.
!  (inout)   CS - The control structure set up by initialize_dyn_unsplit.
!  (inout)   restart_CS - A pointer to the restart control structure.

  type(vardesc) :: vd
  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 would be 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

subroutine initialize_dyn_unsplit(u, v, h, Time, G, GV, param_file, diag, CS, &
                                  restart_CS, Accel_diag, Cont_diag, MIS, &
                                  OBC, update_OBC_CSp, ALE_CSp, setVisc_CSp, &
                                  visc, dirs, ntrunc)
  type(ocean_grid_type),          intent(inout) :: G          !< The ocean's grid structure.
  type(verticalgrid_type),        intent(in)    :: GV         !< The ocean's vertical grid
                                                              !! structure.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                                  intent(inout) :: u          !< The zonal velocity, in m s-1.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                                  intent(inout) :: v          !< The meridional velocity, in m s-1.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) , &
                                  intent(inout) :: h          !< Layer thicknesses, in H
                                                              !! (usually 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(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).

! Arguments: u - The zonal velocity, in m s-1.
!  (inout)   v - The meridional velocity, in m s-1.
!  (inout)   h - The layer thicknesses, in m or kg m-2, depending on whether
!                the Boussinesq approximation is made.
!  (in)      Time - The current model time.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.
!  (in)      diag - A structure that is used to regulate diagnostic output.
!  (inout)   CS - The control structure set up by initialize_dyn_unsplit.
!  (in)      restart_CS - A pointer to the restart control structure.
!  (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.
!  (inout)   Cont_diag - A structure with pointers to various terms in the
!                   continuity equations.
!  (inout)   MIS - The "MOM6 Internal State" structure, used to pass around
!                  pointers to various arrays for diagnostic purposes.
!  (in)      OBC - If open boundary conditions are used, this points to the
!                  ocean_OBC_type that was set up in MOM_initialization.
!  (in)      update_OBC_CSp - If open boundary condition updates are used,
!                  this points to the appropriate control structure.
!  (in)      ALE_CS - This points to the ALE control structure.
!  (in)      setVisc_CSp - This points to the set_visc control structure.
!  (inout)   visc - A structure containing vertical viscosities, bottom drag
!                   viscosities, and related fields.
!  (in)      dirs - A structure containing several relevant directory paths.
!  (in)      ntrunc - A target for the variable that records the number of times
!                     the velocity is truncated (this should be 0).

  !   This subroutine initializes all of the variables that are used by this
  ! dynamic core, including diagnostics and the cpu clocks.
  character(len=40) :: mdl = "MOM_dynamics_unsplit" ! This module's name.
  character(len=48) :: thickness_units, flux_units
  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 get_param(param_file, mdl, "DEBUG", cs%debug, &
                 "If true, write out verbose debugging data.", default=.false.)
  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, param_file, diag, cs%continuity_CSp)
  call coriolisadv_init(time, g, 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, param_file, diag, cs%PressureForce_CSp, &
                          cs%tides_CSp)
  call hor_visc_init(time, g, param_file, diag, cs%hor_visc_CSp)
  call vertvisc_init(mis, time, g, gv, 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)
  cs%id_uh = register_diag_field('ocean_model', 'uh', diag%axesCuL, time, &
      'Zonal Thickness Flux', flux_units, y_cell_method='sum', v_extensive=.true.)
  cs%id_vh = register_diag_field('ocean_model', 'vh', diag%axesCvL, time, &
      'Meridional Thickness Flux', flux_units, x_cell_method='sum', v_extensive=.true.)
  cs%id_CAu = register_diag_field('ocean_model', 'CAu', diag%axesCuL, time, &
      'Zonal Coriolis and Advective Acceleration', 'meter second-2')
  cs%id_CAv = register_diag_field('ocean_model', 'CAv', diag%axesCvL, time, &
      'Meridional Coriolis and Advective Acceleration', 'meter second-2')
  cs%id_PFu = register_diag_field('ocean_model', 'PFu', diag%axesCuL, time, &
      'Zonal Pressure Force Acceleration', 'meter second-2')
  cs%id_PFv = register_diag_field('ocean_model', 'PFv', diag%axesCvL, time, &
      'Meridional Pressure Force Acceleration', 'meter second-2')

  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

subroutine end_dyn_unsplit(CS)
  type(mom_dyn_unsplit_cs), pointer :: CS
!  (inout)    CS - The control structure set up by initialize_dyn_unsplit.

  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