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

integer	id_clock_cor
	CPU time clock IDs. More...

integer	id_clock_pres
	CPU time clock IDs. More...

integer	id_clock_vertvisc
	CPU time clock IDs. More...

integer	id_clock_horvisc
	CPU time clock IDs. More...

integer	id_clock_continuity
	CPU time clock IDs. More...

integer	id_clock_mom_update
	CPU time clock IDs. More...

integer	id_clock_pass
	CPU time clock IDs. More...

integer	id_clock_pass_init
	CPU time clock IDs. More...

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

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

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

Definition at line 510 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).
  
   !   This subroutine initializes all of the variables that are used by this
   ! dynamic core, including diagnostics and the cpu clocks.
  
   ! Local varaibles
   character(len=40) :: mdl = "MOM_dynamics_unsplit_RK2" ! This module's name.
   character(len=48) :: thickness_units, flux_units
   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 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, "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)
   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)
  

References mom_continuity::continuity_init(), mom_coriolisadv::coriolisadv_init(), mom_verticalgrid::get_flux_units(), mom_hor_visc::hor_visc_init(), id_clock_continuity, id_clock_cor, id_clock_horvisc, id_clock_mom_update, id_clock_pass, id_clock_pass_init, id_clock_pres, id_clock_vertvisc, mom_pressureforce::pressureforce_init(), mom_tidal_forcing::tidal_forcing_init(), and mom_vert_friction::vertvisc_init().

Here is the call graph for this function:

◆ 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 465 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.
  

References mom_verticalgrid::get_flux_units().

Here is the call graph for this function:

◆ 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 [Pa].
	p_surf_end	A pointer (perhaps NULL) to the surface pressure at the end of this dynamic step [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 191 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 [Pa].
   real, dimension(:,:),              pointer       :: p_surf_end   !< A pointer (perhaps NULL) to
                                                               !! the surface pressure at the end of
                                                               !! this dynamic step [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].
   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)
   call set_viscous_ml(up, vp, h_av, tv, forces, visc, dt_pred, 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)
  