mom_dynamics_unsplit_rk2 Module Reference

Detailed Description

Time steps the ocean dynamics with an unsplit quasi 2nd order Runge-Kutta scheme.

Data Types
type	mom_dyn_unsplit_rk2_cs
	MOM_dynamics_unsplit_RK2 module control structure. More...

Functions/Subroutines
subroutine, public	step_mom_dyn_unsplit_rk2 (u_in, v_in, h_in, tv, visc, Time_local, dt, forces, p_surf_begin, p_surf_end, uh, vh, uhtr, vhtr, eta_av, G, GV, US, CS, VarMix, MEKE)
	Step the MOM6 dynamics using an unsplit quasi-2nd order Runge-Kutta scheme. More...
subroutine, public	register_restarts_dyn_unsplit_rk2 (HI, GV, param_file, CS, restart_CS)
	Allocate the control structure for this module, allocates memory in it, and registers any auxiliary restart variables that are specific to the unsplit RK2 time stepping scheme. More...
subroutine, public	initialize_dyn_unsplit_rk2 (u, v, h, Time, G, GV, US, param_file, diag, CS, restart_CS, Accel_diag, Cont_diag, MIS, MEKE, OBC, update_OBC_CSp, ALE_CSp, setVisc_CSp, visc, dirs, ntrunc, cont_stencil)
	Initialize parameters and allocate memory associated with the unsplit RK2 dynamics module. More...
subroutine, public	end_dyn_unsplit_rk2 (CS)
	Clean up and deallocate memory associated with the dyn_unsplit_RK2 module. More...

Variables
integer	id_clock_cor
	CPU time clock IDs.
integer	id_clock_pres
	CPU time clock IDs.
integer	id_clock_vertvisc
	CPU time clock IDs.
integer	id_clock_horvisc
	CPU time clock IDs.
integer	id_clock_continuity
	CPU time clock IDs.
integer	id_clock_mom_update
	CPU time clock IDs.
integer	id_clock_pass
	CPU time clock IDs.
integer	id_clock_pass_init
	CPU time clock IDs.

Function/Subroutine Documentation

end_dyn_unsplit_rk2()

subroutine, public mom_dynamics_unsplit_rk2::end_dyn_unsplit_rk2

(

type(mom_dyn_unsplit_rk2_cs), pointer

CS

)

Clean up and deallocate memory associated with the dyn_unsplit_RK2 module.

Parameters

cs	dyn_unsplit_RK2 control structure that will be deallocated in this subroutine.

Definition at line 666 of file MOM_dynamics_unsplit_RK2.F90.

  type(MOM_dyn_unsplit_RK2_CS), pointer :: CS !< dyn_unsplit_RK2 control structure that
                                              !! will be deallocated in this subroutine.

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

  deallocate(cs)

initialize_dyn_unsplit_rk2()

subroutine, public mom_dynamics_unsplit_rk2::initialize_dyn_unsplit_rk2	(	real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(inout)	u,
		real, dimension( g %isd: g %ied, g %jsdb: g %jedb, g %ke), intent(inout)	v,
		real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(inout)	h,
		type(time_type), intent(in), target	Time,
		type(ocean_grid_type), intent(inout)	G,
		type(verticalgrid_type), intent(in)	GV,
		type(unit_scale_type), intent(in)	US,
		type(param_file_type), intent(in)	param_file,
		type(diag_ctrl), intent(inout), target	diag,
		type(mom_dyn_unsplit_rk2_cs), pointer	CS,
		type(mom_restart_cs), pointer	restart_CS,
		type(accel_diag_ptrs), intent(inout), target	Accel_diag,
		type(cont_diag_ptrs), intent(inout), target	Cont_diag,
		type(ocean_internal_state), intent(inout)	MIS,
		type(meke_type), pointer	MEKE,
		type(ocean_obc_type), pointer	OBC,
		type(update_obc_cs), pointer	update_OBC_CSp,
		type(ale_cs), pointer	ALE_CSp,
		type(set_visc_cs), pointer	setVisc_CSp,
		type(vertvisc_type), intent(inout)	visc,
		type(directories), intent(in)	dirs,
		integer, intent(inout), target	ntrunc,
		integer, intent(out), optional	cont_stencil
	)

Initialize parameters and allocate memory associated with the unsplit RK2 dynamics module.

Parameters

[in,out]	g	The ocean's grid structure.
[in]	gv	The ocean's vertical grid structure.
[in]	us	A dimensional unit scaling type
[in,out]	u	The zonal velocity [L T-1 ~> m s-1].
[in,out]	v	The meridional velocity [L T-1 ~> m s-1].
[in,out]	h	Layer thicknesses [H ~> m or kg m-2]
[in]	time	The current model time.
[in]	param_file	A structure to parse for run-time parameters.
[in,out]	diag	A structure that is used to regulate diagnostic output.
	cs	The control structure set up by initialize_dyn_unsplit_RK2.
	restart_cs	A pointer to the restart control structure.
[in,out]	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.
[in,out]	cont_diag	A structure with pointers to various terms in the continuity equations.
[in,out]	mis	The "MOM6 Internal State" structure, used to pass around pointers to various arrays for diagnostic purposes.
	meke	MEKE data
	obc	If open boundary conditions are used, this points to the ocean_OBC_type that was set up in MOM_initialization.
	update_obc_csp	If open boundary condition updates are used, this points to the appropriate control structure.
	ale_csp	This points to the ALE control structure.
	setvisc_csp	This points to the set_visc control structure.
[in,out]	visc	A structure containing vertical viscosities, bottom drag viscosities, and related fields.
[in]	dirs	A structure containing several relevant directory paths.
[in,out]	ntrunc	A target for the variable that records the number of times the velocity is truncated (this should be 0).
[out]	cont_stencil	The stencil for thickness from the continuity solver.

Definition at line 508 of file MOM_dynamics_unsplit_RK2.F90.

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

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

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

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

  cs%diag => diag

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

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

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

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

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

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

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

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

register_restarts_dyn_unsplit_rk2()

subroutine, public mom_dynamics_unsplit_rk2::register_restarts_dyn_unsplit_rk2	(	type(hor_index_type), intent(in)	HI,
		type(verticalgrid_type), intent(in)	GV,
		type(param_file_type), intent(in)	param_file,
		type(mom_dyn_unsplit_rk2_cs), pointer	CS,
		type(mom_restart_cs), pointer	restart_CS
	)

Allocate the control structure for this module, allocates memory in it, and registers any auxiliary restart variables that are specific to the unsplit RK2 time stepping scheme.

All variables registered here should have the ability to be recreated if they are not present in a restart file.

Parameters

[in]	hi	A horizontal index type structure.
[in]	gv	The ocean's vertical grid structure.
[in]	param_file	A structure to parse for run-time parameters.
	cs	The control structure set up by initialize_dyn_unsplit_RK2.
	restart_cs	A pointer to the restart control structure.

Definition at line 463 of file MOM_dynamics_unsplit_RK2.F90.

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

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

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

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

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

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

step_mom_dyn_unsplit_rk2()

subroutine, public mom_dynamics_unsplit_rk2::step_mom_dyn_unsplit_rk2	(	real, dimension(szib_(g),szj_(g),szk_(g)), intent(inout)	u_in,
		real, dimension(szi_(g),szjb_(g),szk_(g)), intent(inout)	v_in,
		real, dimension(szi_(g),szj_(g),szk_(g)), intent(inout)	h_in,
		type(thermo_var_ptrs), intent(in)	tv,
		type(vertvisc_type), intent(inout)	visc,
		type(time_type), intent(in)	Time_local,
		real, intent(in)	dt,
		type(mech_forcing), intent(in)	forces,
		real, dimension(:,:), pointer	p_surf_begin,
		real, dimension(:,:), pointer	p_surf_end,
		real, dimension(szib_(g),szj_(g),szk_(g)), intent(inout)	uh,
		real, dimension(szi_(g),szjb_(g),szk_(g)), intent(inout)	vh,
		real, dimension(szib_(g),szj_(g),szk_(g)), intent(inout)	uhtr,
		real, dimension(szi_(g),szjb_(g),szk_(g)), intent(inout)	vhtr,
		real, dimension(szi_(g),szj_(g)), intent(out)	eta_av,
		type(ocean_grid_type), intent(inout)	G,
		type(verticalgrid_type), intent(in)	GV,
		type(unit_scale_type), intent(in)	US,
		type(mom_dyn_unsplit_rk2_cs), pointer	CS,
		type(varmix_cs), pointer	VarMix,
		type(meke_type), pointer	MEKE
	)

Step the MOM6 dynamics using an unsplit quasi-2nd order Runge-Kutta scheme.

Parameters

[in,out]	g	The ocean's grid structure.
[in]	gv	The ocean's vertical grid structure.
[in]	us	A dimensional unit scaling type
[in,out]	u_in	The input and output zonal velocity [L T-1 ~> m s-1].
[in,out]	v_in	The input and output meridional velocity [L T-1 ~> m s-1].
[in,out]	h_in	The input and output layer thicknesses, [H ~> m or kg m-2], depending on whether the Boussinesq approximation is made.
[in]	tv	A structure pointing to various thermodynamic variables.
[in,out]	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 baroclinic dynamics time step [T ~> s].
[in]	forces	A structure with the driving mechanical forces
	p_surf_begin	A pointer (perhaps NULL) to the surface pressure at the beginning of this dynamic step [R L2 T-2 ~> Pa].
	p_surf_end	A pointer (perhaps NULL) to the surface pressure at the end of this dynamic step [R L2 T-2 ~> Pa].
[in,out]	uh	The zonal volume or mass transport [H L2 T-1 ~> m3 s-1 or kg s-1].
[in,out]	vh	The meridional volume or mass transport [H L2 T-1 ~> m3 s-1 or kg s-1].
[in,out]	uhtr	The accumulated zonal volume or mass transport since the last tracer advection [H L2 ~> m3 or kg].
[in,out]	vhtr	The accumulated meridional volume or mass transport since the last tracer advection [H L2 ~> m3 or kg].
[out]	eta_av	The time-mean free surface height or column mass [H ~> m or kg m-2].
	cs	The control structure set up by initialize_dyn_unsplit_RK2.
	varmix	A pointer to a structure with fields that specify the spatially variable viscosities.
	meke	A pointer to a structure containing fields related to the Mesoscale Eddy Kinetic Energy.

Definition at line 194 of file MOM_dynamics_unsplit_RK2.F90.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  if (dyn_p_surf) deallocate(p_surf)

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

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