MOM_dynamics_unsplit_RK2.F90

Go to the documentation of this file.

module mom_dynamics_unsplit_rk2

!***********************************************************************
!*                   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 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 : 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
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, 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_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_rk2_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.

  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 :: 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(tidal_forcing_cs), pointer :: tides_csp => null()
  type(update_obc_cs), pointer :: update_obc_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_rk2_cs


public step_mom_dyn_unsplit_rk2, register_restarts_dyn_unsplit_rk2
public initialize_dyn_unsplit_rk2, end_dyn_unsplit_rk2

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

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

subroutine step_mom_dyn_unsplit_rk2(u_in, v_in, h_in, 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_in    !< The input and output zonal
                                                              !! velocity, in m s-1.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                                     intent(inout) :: v_in    !< The input and output meridional
                                                              !! velocity, in m s-1.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  &
                                     intent(inout) :: h_in    !< The input and output layer
                                              !! thicknesses, in 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,
                                                              !! 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    !< The 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_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.
! Arguments: u_in - The input and output zonal velocity, in m s-1.
!  (inout)   v_in - The input and output meridional velocity, in m s-1.
!  (inout)   h_in - 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_RK2.
!  (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
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: vp
  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 * 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 => fluxes%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)
  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_in, v_in, h_in, cs%diffu, cs%diffv, meke, varmix, &
                            g, gv, cs%hor_visc_CSp)
  call cpu_clock_end(id_clock_horvisc)
  call disable_averaging(cs%diag)
  call cpu_clock_begin(id_clock_pass)
  call pass_vector(cs%diffu, cs%diffv, g%Domain)
  call cpu_clock_end(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, 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)

! 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, 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, &
                     cs%PressureForce_CSp, cs%ALE_CSp, p_surf)
  call cpu_clock_end(id_clock_pres)
  call cpu_clock_begin(id_clock_pass)
  call pass_vector(cs%PFu, cs%PFv, g%Domain)
  call pass_vector(cs%CAu, cs%CAv, g%Domain)
  call cpu_clock_end(id_clock_pass)

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

 ! 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_averaging(dt, time_local, cs%diag)
  call set_viscous_ml(up, vp, h_av, tv, fluxes, visc, dt_pred, g, gv, &
                      cs%set_visc_CSp)
  call disable_averaging(cs%diag)
  call vertvisc_coef(up, vp, h_av, fluxes, visc, dt_pred, g, gv, &
                     cs%vertvisc_CSp, cs%OBC)
  call vertvisc(up, vp, h_av, fluxes, visc, dt_pred, 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[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, 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)

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

! 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, 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_averaging(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, fluxes, visc, dt, g, gv, &
                     cs%vertvisc_CSp, cs%OBC)
  call vertvisc(up, vp, h_av, fluxes, visc, dt, cs%OBC, cs%ADp, cs%CDp, &
                g, gv, cs%vertvisc_CSp, cs%taux_bot, cs%tauy_bot)
  call vertvisc_coef(u_in, v_in, h_av, fluxes, visc, dt, g, gv, &
                     cs%vertvisc_CSp, cs%OBC)
  call vertvisc(u_in, v_in, 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(up, vp, g%Domain)
  call pass_vector(u_in, v_in, g%Domain)
  call cpu_clock_end(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, cs%continuity_CSp, obc=cs%OBC)
  call cpu_clock_end(id_clock_continuity)
  call cpu_clock_begin(id_clock_pass)
  call pass_var(h_in, g%Domain)
  call pass_vector(uh, vh, g%Domain)
  call cpu_clock_end(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)
    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)

!   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

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

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.

! 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_RK2.
!  (inout)   restart_CS - A pointer to the restart control structure.

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

subroutine initialize_dyn_unsplit_rk2(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_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(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_RK2.
!  (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_RK2" ! 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_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 get_param(param_file, mdl, "BE", cs%be, &
                 "If SPLIT is true, BE determines the relative weighting \n"//&
                 "of a  2nd-order Runga-Kutta baroclinic time stepping \n"//&
                 "scheme (0.5) and a backward Euler scheme (1) that is \n"//&
                 "used for the Coriolis and inertial terms.  BE may be \n"//&
                 "from 0.5 to 1, but instability may occur near 0.5. \n"//&
                 "BE is also applicable if SPLIT is false and USE_RK2 \n"//&
                 "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 \n"//&
                 "controls the extent to which the treatment of gravity \n"//&
                 "waves is forward-backward (0) or simulated backward \n"//&
                 "Euler (1).  0 is almost always used.\n"//&
                 "If SPLIT is false and USE_RK2 is true, BEGW can be \n"//&
                 "between 0 and 0.5 to damp gravity waves.", &
                 units="nondim", default=0.0)
  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_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)
  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_rk2

subroutine end_dyn_unsplit_rk2(CS)
  type(mom_dyn_unsplit_rk2_cs), pointer :: CS
!  (inout)   CS - The control structure set up by initialize_dyn_unsplit_RK2.

  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