mom_dynamics_unsplit Module Reference

Detailed Description

Time steps the ocean dynamics with an unsplit quasi 3rd order scheme.

Data Types
type	mom_dyn_unsplit_cs
	MOM_dynamics_unsplit module control structure. More...

Functions/Subroutines
subroutine, public	step_mom_dyn_unsplit (u, v, h, tv, visc, Time_local, dt, forces, p_surf_begin, p_surf_end, uh, vh, uhtr, vhtr, eta_av, G, GV, US, CS, VarMix, MEKE, Waves)
	Step the MOM6 dynamics using an unsplit mixed 2nd order (for continuity) and 3rd order (for the inviscid momentum equations) order scheme. More...
subroutine, public	register_restarts_dyn_unsplit (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 time stepping scheme. More...
subroutine, public	initialize_dyn_unsplit (u, v, h, Time, G, GV, US, param_file, diag, CS, restart_CS, Accel_diag, Cont_diag, MIS, MEKE, OBC, update_OBC_CSp, ALE_CSp, setVisc_CSp, visc, dirs, ntrunc, cont_stencil)
	Initialize parameters and allocate memory associated with the unsplit dynamics module. More...
subroutine, public	end_dyn_unsplit (CS)
	Clean up and deallocate memory associated with the unsplit dynamics 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_continuity
	CPU time clock IDs.
integer	id_clock_horvisc
	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()

subroutine, public mom_dynamics_unsplit::end_dyn_unsplit

(

type(mom_dyn_unsplit_cs), pointer

CS

)

Clean up and deallocate memory associated with the unsplit dynamics module.

Parameters

cs	unsplit dynamics control structure that will be deallocated in this subroutine.

Definition at line 708 of file MOM_dynamics_unsplit.F90.

  type(MOM_dyn_unsplit_CS), pointer :: CS !< unsplit dynamics control structure that
                                          !! will be deallocated in this subroutine.
 
  dealloc_(cs%diffu) ; dealloc_(cs%diffv)
  dealloc_(cs%CAu)   ; dealloc_(cs%CAv)
  dealloc_(cs%PFu)   ; dealloc_(cs%PFv)
 
  deallocate(cs)

initialize_dyn_unsplit()

subroutine, public mom_dynamics_unsplit::initialize_dyn_unsplit	(	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_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 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.
	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 556 of file MOM_dynamics_unsplit.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_CS),       pointer       :: CS         !< The control structure set up
                                                              !! by initialize_dyn_unsplit.
  type(MOM_restart_CS),           pointer       :: restart_CS !< A pointer to the restart control
                                                              !!structure.
  type(accel_diag_ptrs),  target, intent(inout) :: Accel_diag !< A set of pointers to the various
                                     !! accelerations in the momentum equations, which can be used
                                     !! for later derived diagnostics, like energy budgets.
  type(cont_diag_ptrs),   target, intent(inout) :: Cont_diag  !< A structure with pointers to
                                                              !! various terms in the continuity
                                                              !! equations.
  type(ocean_internal_state),     intent(inout) :: MIS        !< The "MOM6 Internal State"
                                                   !! structure, used to pass around pointers
                                                   !! to various arrays for diagnostic purposes.
  type(MEKE_type),                pointer       :: MEKE !< MEKE data
  type(ocean_OBC_type),           pointer       :: OBC        !< If open boundary conditions are
                                                       !! used, this points to the ocean_OBC_type
                                                       !! that was set up in MOM_initialization.
  type(update_OBC_CS),            pointer       :: update_OBC_CSp !< If open boundary condition
                                                            !! updates are used, this points to
                                                            !! the appropriate control structure.
  type(ALE_CS),                   pointer       :: ALE_CSp    !< This points to the ALE control
                                                              !! structure.
  type(set_visc_CS),              pointer       :: setVisc_CSp !< This points to the set_visc
                                                               !! control structure.
  type(vertvisc_type),            intent(inout) :: visc       !< A structure containing vertical
                                                              !! viscosities, bottom drag
                                                              !! viscosities, and related fields.
  type(directories),              intent(in)    :: dirs       !< A structure containing several
                                                              !! relevant directory paths.
  integer, target,                intent(inout) :: ntrunc     !< A target for the variable that
                                                        !! records the number of times the velocity
                                                        !! is truncated (this should be 0).
  integer,              optional, intent(out)   :: cont_stencil !< The stencil for thickness
                                                                !! from the continuity solver.
 
  !   This subroutine initializes all of the variables that are used by this
  ! dynamic core, including diagnostics and the cpu clocks.
 
  ! Local variables
  character(len=40) :: mdl = "MOM_dynamics_unsplit" ! This module's name.
  character(len=48) :: thickness_units, flux_units
  ! This include declares and sets the variable "version".
# include "version_variable.h"
  real :: H_convert
  logical :: use_tides
  integer :: isd, ied, jsd, jed, nz, IsdB, IedB, JsdB, JedB
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = g%ke
  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
 
  if (.not.associated(cs)) call mom_error(fatal, &
      "initialize_dyn_unsplit called with an unassociated control structure.")
  if (cs%module_is_initialized) then
    call mom_error(warning, "initialize_dyn_unsplit called with a control "// &
                            "structure that has already been initialized.")
    return
  endif
  cs%module_is_initialized = .true.
 
  cs%diag => diag
 
  call log_version(param_file, mdl, version, "")
  call get_param(param_file, mdl, "FIX_UNSPLIT_DT_VISC_BUG", cs%use_correct_dt_visc, &
                 "If true, use the correct timestep in the viscous terms applied in the first "//&
                 "predictor step with the unsplit time stepping scheme, and in the calculation "//&
                 "of the turbulent mixed layer properties for viscosity with unsplit or "//&
                 "unsplit_RK2.", default=.true.)
  call get_param(param_file, mdl, "DEBUG", cs%debug, &
                 "If true, write out verbose debugging data.", &
                 default=.false., debuggingparam=.true.)
  call get_param(param_file, mdl, "TIDES", use_tides, &
                 "If true, apply tidal momentum forcing.", default=.false.)
 
  allocate(cs%taux_bot(isdb:iedb,jsd:jed)) ; cs%taux_bot(:,:) = 0.0
  allocate(cs%tauy_bot(isd:ied,jsdb:jedb)) ; cs%tauy_bot(:,:) = 0.0
 
  mis%diffu => cs%diffu ; mis%diffv => cs%diffv
  mis%PFu => cs%PFu ; mis%PFv => cs%PFv
  mis%CAu => cs%CAu ; mis%CAv => cs%CAv
 
  cs%ADp => accel_diag ; cs%CDp => cont_diag
  accel_diag%diffu => cs%diffu ; accel_diag%diffv => cs%diffv
  accel_diag%PFu => cs%PFu ; accel_diag%PFv => cs%PFv
  accel_diag%CAu => cs%CAu ; accel_diag%CAv => cs%CAv
 
  call continuity_init(time, g, gv, us, param_file, diag, cs%continuity_CSp)
  if (present(cont_stencil)) cont_stencil = continuity_stencil(cs%continuity_CSp)
  call coriolisadv_init(time, g, gv, us, param_file, diag, cs%ADp, cs%CoriolisAdv_CSp)
  if (use_tides) call tidal_forcing_init(time, g, param_file, cs%tides_CSp)
  call pressureforce_init(time, g, gv, us, param_file, diag, cs%PressureForce_CSp, &
                          cs%tides_CSp)
  call hor_visc_init(time, g, us, param_file, diag, cs%hor_visc_CSp, meke)
  call vertvisc_init(mis, time, g, gv, us, param_file, diag, cs%ADp, dirs, &
                     ntrunc, cs%vertvisc_CSp)
  if (.not.associated(setvisc_csp)) call mom_error(fatal, &
    "initialize_dyn_unsplit called with setVisc_CSp unassociated.")
  cs%set_visc_CSp => setvisc_csp
 
  if (associated(ale_csp)) cs%ALE_CSp => ale_csp
  if (associated(obc)) cs%OBC => obc
  if (associated(update_obc_csp)) cs%update_OBC_CSp => update_obc_csp
 
  flux_units = get_flux_units(gv)
  h_convert = gv%H_to_m ; if (.not.gv%Boussinesq) h_convert = gv%H_to_kg_m2
  cs%id_uh = register_diag_field('ocean_model', 'uh', diag%axesCuL, time, &
      'Zonal Thickness Flux', flux_units, y_cell_method='sum', v_extensive=.true., &
      conversion=h_convert*us%L_to_m**2*us%s_to_T)
  cs%id_vh = register_diag_field('ocean_model', 'vh', diag%axesCvL, time, &
      'Meridional Thickness Flux', flux_units, x_cell_method='sum', v_extensive=.true., &
      conversion=h_convert*us%L_to_m**2*us%s_to_T)
  cs%id_CAu = register_diag_field('ocean_model', 'CAu', diag%axesCuL, time, &
      'Zonal Coriolis and Advective Acceleration', 'm s-2', &
      conversion=us%L_T2_to_m_s2)
  cs%id_CAv = register_diag_field('ocean_model', 'CAv', diag%axesCvL, time, &
      'Meridional Coriolis and Advective Acceleration', 'm s-2', &
      conversion=us%L_T2_to_m_s2)
  cs%id_PFu = register_diag_field('ocean_model', 'PFu', diag%axesCuL, time, &
      'Zonal Pressure Force Acceleration', 'm s-2', &
      conversion=us%L_T2_to_m_s2)
  cs%id_PFv = register_diag_field('ocean_model', 'PFv', diag%axesCvL, time, &
      'Meridional Pressure Force Acceleration', 'm s-2', &
      conversion=us%L_T2_to_m_s2)
 
  id_clock_cor = cpu_clock_id('(Ocean Coriolis & mom advection)', grain=clock_module)
  id_clock_continuity = cpu_clock_id('(Ocean continuity equation)', grain=clock_module)
  id_clock_pres = cpu_clock_id('(Ocean pressure force)', grain=clock_module)
  id_clock_vertvisc = cpu_clock_id('(Ocean vertical viscosity)', grain=clock_module)
  id_clock_horvisc = cpu_clock_id('(Ocean horizontal viscosity)', grain=clock_module)
  id_clock_mom_update = cpu_clock_id('(Ocean momentum increments)', grain=clock_module)
  id_clock_pass = cpu_clock_id('(Ocean message passing)', grain=clock_module)
  id_clock_pass_init = cpu_clock_id('(Ocean init message passing)', grain=clock_routine)
 

register_restarts_dyn_unsplit()

subroutine, public mom_dynamics_unsplit::register_restarts_dyn_unsplit	(	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_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 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.
	restart_cs	A pointer to the restart control structure.

Definition at line 517 of file MOM_dynamics_unsplit.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_CS),  pointer    :: CS         !< The control structure set up by
                                                      !! initialize_dyn_unsplit.
  type(MOM_restart_CS),      pointer    :: restart_CS !< A pointer to the restart control structure.
 
  ! Local arguments
  character(len=40)  :: mdl = "MOM_dynamics_unsplit" ! This module's name.
  character(len=48) :: thickness_units, flux_units
  integer :: isd, ied, jsd, jed, nz, IsdB, IedB, JsdB, JedB
  isd = hi%isd ; ied = hi%ied ; jsd = hi%jsd ; jed = hi%jed ; nz = gv%ke
  isdb = hi%IsdB ; iedb = hi%IedB ; jsdb = hi%JsdB ; jedb = hi%JedB
 
  ! This is where a control structure that is specific to this module is allocated.
  if (associated(cs)) then
    call mom_error(warning, "register_restarts_dyn_unsplit called with an associated "// &
                             "control structure.")
    return
  endif
  allocate(cs)
 
  alloc_(cs%diffu(isdb:iedb,jsd:jed,nz)) ; cs%diffu(:,:,:) = 0.0
  alloc_(cs%diffv(isd:ied,jsdb:jedb,nz)) ; cs%diffv(:,:,:) = 0.0
  alloc_(cs%CAu(isdb:iedb,jsd:jed,nz)) ; cs%CAu(:,:,:) = 0.0
  alloc_(cs%CAv(isd:ied,jsdb:jedb,nz)) ; cs%CAv(:,:,:) = 0.0
  alloc_(cs%PFu(isdb:iedb,jsd:jed,nz)) ; cs%PFu(:,:,:) = 0.0
  alloc_(cs%PFv(isd:ied,jsdb:jedb,nz)) ; cs%PFv(:,:,:) = 0.0
 
  thickness_units = get_thickness_units(gv)
  flux_units = get_flux_units(gv)
 
!  No extra restart fields are needed with this time stepping scheme.
 

step_mom_dyn_unsplit()

subroutine, public mom_dynamics_unsplit::step_mom_dyn_unsplit	(	real, dimension(szib_(g),szj_(g),szk_(g)), intent(inout)	u,
		real, dimension(szi_(g),szjb_(g),szk_(g)), intent(inout)	v,
		real, dimension(szi_(g),szj_(g),szk_(g)), intent(inout)	h,
		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_cs), pointer	CS,
		type(varmix_cs), pointer	VarMix,
		type(meke_type), pointer	MEKE,
		type(wave_parameters_cs), optional, pointer	Waves
	)

Step the MOM6 dynamics using an unsplit mixed 2nd order (for continuity) and 3rd order (for the inviscid momentum equations) order 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	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]	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 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 start 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.
	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.
	waves	A pointer to a structure containing fields related to the surface wave conditions

Definition at line 189 of file MOM_dynamics_unsplit.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(thermo_var_ptrs),   intent(in)    :: tv     !< A structure pointing to various
                                                   !! thermodynamic variables.
  type(vertvisc_type),     intent(inout) :: visc   !< A structure containing vertical
                                 !! viscosities, bottom drag viscosities, and related fields.
  type(time_type),         intent(in)    :: Time_local   !< The model time at the end
                                                         !! of the time step.
  real,                    intent(in)    :: dt     !< The dynamics time step [T ~> s].
  type(mech_forcing),      intent(in)    :: forces !< A structure with the driving mechanical forces
  real, dimension(:,:),    pointer       :: p_surf_begin !< A pointer (perhaps NULL) to the surface
                                                   !! pressure at the start of this dynamic step [R L2 T-2 ~> Pa].
  real, dimension(:,:),    pointer       :: p_surf_end   !< A pointer (perhaps NULL) to the surface
                                                   !! pressure at the end of this dynamic step [R L2 T-2 ~> Pa].
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: uh !< The zonal volume or mass transport
                                                   !! [H L2 T-1 ~> m3 s-1 or kg s-1].
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: vh !< The meridional volume or mass
                                                   !! transport [H L2 T-1 ~> m3 s-1 or kg s-1].
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: uhtr !< The accumulated zonal volume or mass
                                                   !! transport since the last tracer advection [H L2 ~> m3 or kg].
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: vhtr !< The accumulated meridional volume or mass
                                                   !! transport since the last tracer advection [H L2 ~> m3 or kg].
  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: eta_av !< The time-mean free surface height or
                                                   !! column mass [H ~> m or kg m-2].
  type(MOM_dyn_unsplit_CS), pointer      :: CS     !< The control structure set up by
                                                   !! initialize_dyn_unsplit.
  type(VarMix_CS),         pointer       :: VarMix !< A pointer to a structure with fields
                                 !! that specify the spatially variable viscosities.
  type(MEKE_type),         pointer       :: MEKE   !< A pointer to a structure containing
                                 !! fields related to the Mesoscale Eddy Kinetic Energy.
  type(wave_parameters_CS), optional, pointer :: Waves !< A pointer to a structure containing
                                 !! fields related to the surface wave conditions
 
  ! Local variables
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_av, hp ! Prediced or averaged layer thicknesses [H ~> m or kg m-2]
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: up, upp ! Predicted zonal velocities [L T-1 ~> m s-1]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: vp, vpp ! Predicted meridional velocities [L T-1 ~> m s-1]
  real, dimension(:,:), pointer :: p_surf => null()
  real :: dt_pred   ! The time step for the predictor part of the baroclinic time stepping [T ~> s].
  real :: dt_visc   ! The time step for a part of the update due to viscosity [T ~> s].
  logical :: dyn_p_surf
  integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
  dt_pred = dt / 3.0
 
  h_av(:,:,:) = 0; hp(:,:,:) = 0
  up(:,:,:) = 0; upp(:,:,:) = 0
  vp(:,:,:) = 0; vpp(:,:,:) = 0
 
  dyn_p_surf = associated(p_surf_begin) .and. associated(p_surf_end)
  if (dyn_p_surf) then
    call safe_alloc_ptr(p_surf,g%isd,g%ied,g%jsd,g%jed) ; p_surf(:,:) = 0.0
  else
    p_surf => forces%p_surf
  endif
 
! Matsuno's third order accurate three step scheme is used to step
! all of the fields except h.  h is stepped separately.
 
  if (cs%debug) then
    call mom_state_chksum("Start First Predictor ", u, v, h, uh, vh, g, gv, us)
  endif
 
! diffu = horizontal viscosity terms (u,h)
  call enable_averages(dt, time_local, cs%diag)
  call cpu_clock_begin(id_clock_horvisc)
  call horizontal_viscosity(u, v, h, cs%diffu, cs%diffv, meke, varmix, g, gv, us, cs%hor_visc_CSp)
  call cpu_clock_end(id_clock_horvisc)
  call disable_averaging(cs%diag)
 
! uh = u*h
! hp = h + dt/2 div . uh
  call cpu_clock_begin(id_clock_continuity)
  call continuity(u, v, h, hp, uh, vh, dt*0.5, g, gv, us, cs%continuity_CSp, obc=cs%OBC)
  call cpu_clock_end(id_clock_continuity)
  call pass_var(hp, g%Domain, clock=id_clock_pass)
  call pass_vector(uh, vh, g%Domain, clock=id_clock_pass)
 
  call enable_averages(0.5*dt, time_local-real_to_time(0.5*us%T_to_s*dt), cs%diag)
!   Here the first half of the thickness fluxes are offered for averaging.
  if (cs%id_uh > 0) call post_data(cs%id_uh, uh, cs%diag)
  if (cs%id_vh > 0) call post_data(cs%id_vh, vh, cs%diag)
  call disable_averaging(cs%diag)
 
! h_av = (h + hp)/2
! u = u + dt diffu
  call cpu_clock_begin(id_clock_mom_update)
  do k=1,nz
    do j=js-2,je+2 ; do i=is-2,ie+2
      h_av(i,j,k) = (h(i,j,k) + hp(i,j,k)) * 0.5
    enddo ; enddo
    do j=js,je ; do i=isq,ieq
      u(i,j,k) = u(i,j,k) + dt * cs%diffu(i,j,k) * g%mask2dCu(i,j)
    enddo ; enddo
    do j=jsq,jeq ; do i=is,ie
      v(i,j,k) = v(i,j,k) + dt * cs%diffv(i,j,k) * g%mask2dCv(i,j)
    enddo ; enddo
    do j=js-2,je+2 ; do i=isq-2,ieq+2
      uhtr(i,j,k) = uhtr(i,j,k) + 0.5*dt*uh(i,j,k)
    enddo ; enddo
    do j=jsq-2,jeq+2 ; do i=is-2,ie+2
      vhtr(i,j,k) = vhtr(i,j,k) + 0.5*dt*vh(i,j,k)
    enddo ; enddo
  enddo
  call cpu_clock_end(id_clock_mom_update)
  call pass_vector(u, v, g%Domain, clock=id_clock_pass)
 
! CAu = -(f+zeta)/h_av vh + d/dx KE
  call cpu_clock_begin(id_clock_cor)
  call coradcalc(u, v, h_av, uh, vh, cs%CAu, cs%CAv, cs%OBC, cs%ADp, &
                 g, gv, us, cs%CoriolisAdv_CSp)
  call cpu_clock_end(id_clock_cor)
 
! PFu = d/dx M(h_av,T,S)
  call cpu_clock_begin(id_clock_pres)
  if (dyn_p_surf) then ; do j=js-2,je+2 ; do i=is-2,ie+2
    p_surf(i,j) = 0.75*p_surf_begin(i,j) + 0.25*p_surf_end(i,j)
  enddo ; enddo ; endif
  call pressureforce(h_av, tv, cs%PFu, cs%PFv, g, gv, us, &
                     cs%PressureForce_CSp, cs%ALE_CSp, p_surf)
  call cpu_clock_end(id_clock_pres)
 
  if (associated(cs%OBC)) then; if (cs%OBC%update_OBC) then
    call update_obc_data(cs%OBC, g, gv, us, tv, h, cs%update_OBC_CSp, time_local)
  endif; endif
  if (associated(cs%OBC)) then
    call open_boundary_zero_normal_flow(cs%OBC, g, cs%PFu, cs%PFv)
    call open_boundary_zero_normal_flow(cs%OBC, g, cs%CAu, cs%CAv)
  endif
 
! up = u + dt_pred * (PFu + CAu)
  call cpu_clock_begin(id_clock_mom_update)
  do k=1,nz ; do j=js,je ; do i=isq,ieq
    up(i,j,k) = g%mask2dCu(i,j) * (u(i,j,k) + dt_pred * (cs%PFu(i,j,k) + cs%CAu(i,j,k)))
  enddo ; enddo ; enddo
  do k=1,nz ; do j=jsq,jeq ; do i=is,ie
    vp(i,j,k) = g%mask2dCv(i,j) * (v(i,j,k) + dt_pred * (cs%PFv(i,j,k) + cs%CAv(i,j,k)))
  enddo ; enddo ; enddo
  call cpu_clock_end(id_clock_mom_update)
 
  if (cs%debug) then
    call mom_state_chksum("Predictor 1", up, vp, h_av, uh, vh, g, gv, us)
    call mom_accel_chksum("Predictor 1 accel", cs%CAu, cs%CAv, cs%PFu, cs%PFv, &
                          cs%diffu, cs%diffv, g, gv, us)
  endif
 
 ! up <- up + dt/2 d/dz visc d/dz up
  call cpu_clock_begin(id_clock_vertvisc)
  call enable_averages(dt, time_local, cs%diag)
  dt_visc = 0.5*dt ; if (cs%use_correct_dt_visc) dt_visc = dt
  call set_viscous_ml(u, v, h_av, tv, forces, visc, dt_visc, g, gv, us, cs%set_visc_CSp)
  call disable_averaging(cs%diag)
 
  dt_visc = 0.5*dt ; if (cs%use_correct_dt_visc) dt_visc = dt_pred
  call vertvisc_coef(up, vp, h_av, forces, visc, dt_visc, g, gv, us, cs%vertvisc_CSp, cs%OBC)
  call vertvisc(up, vp, h_av, forces, visc, dt_visc, cs%OBC, cs%ADp, cs%CDp, &
                g, gv, us, cs%vertvisc_CSp, waves=waves)
  call cpu_clock_end(id_clock_vertvisc)
  call pass_vector(up, vp, g%Domain, clock=id_clock_pass)
 
! uh = up * hp
! h_av = hp + dt/2 div . uh
  call cpu_clock_begin(id_clock_continuity)
  call continuity(up, vp, hp, h_av, uh, vh, (0.5*dt), g, gv, us, cs%continuity_CSp, obc=cs%OBC)
  call cpu_clock_end(id_clock_continuity)
  call pass_var(h_av, g%Domain, clock=id_clock_pass)
  call pass_vector(uh, vh, g%Domain, clock=id_clock_pass)
 
! h_av <- (hp + h_av)/2
  do k=1,nz ; do j=js-2,je+2 ; do i=is-2,ie+2
    h_av(i,j,k) = (hp(i,j,k) + h_av(i,j,k)) * 0.5
  enddo ; enddo ; enddo
 
! CAu = -(f+zeta(up))/h_av vh + d/dx KE(up)
  call cpu_clock_begin(id_clock_cor)
  call coradcalc(up, vp, h_av, uh, vh, cs%CAu, cs%CAv, cs%OBC, cs%ADp, &
                 g, gv, us, cs%CoriolisAdv_CSp)
  call cpu_clock_end(id_clock_cor)
 
! PFu = d/dx M(h_av,T,S)
  call cpu_clock_begin(id_clock_pres)
  if (dyn_p_surf) then ; do j=js-2,je+2 ; do i=is-2,ie+2
    p_surf(i,j) = 0.25*p_surf_begin(i,j) + 0.75*p_surf_end(i,j)
  enddo ; enddo ; endif
  call pressureforce(h_av, tv, cs%PFu, cs%PFv, g, gv, us, &
                     cs%PressureForce_CSp, cs%ALE_CSp, p_surf)
  call cpu_clock_end(id_clock_pres)
 
  if (associated(cs%OBC)) then; if (cs%OBC%update_OBC) then
    call update_obc_data(cs%OBC, g, gv, us, tv, h, cs%update_OBC_CSp, time_local)
  endif; endif
  if (associated(cs%OBC)) then
    call open_boundary_zero_normal_flow(cs%OBC, g, cs%PFu, cs%PFv)
    call open_boundary_zero_normal_flow(cs%OBC, g, cs%CAu, cs%CAv)
  endif
 
! upp = u + dt/2 * ( PFu + CAu )
  call cpu_clock_begin(id_clock_mom_update)
  do k=1,nz ; do j=js,je ; do i=isq,ieq
    upp(i,j,k) = g%mask2dCu(i,j) * (u(i,j,k) + dt * 0.5 * (cs%PFu(i,j,k) + cs%CAu(i,j,k)))
  enddo ; enddo ; enddo
  do k=1,nz ; do j=jsq,jeq ; do i=is,ie
    vpp(i,j,k) = g%mask2dCv(i,j) * (v(i,j,k) + dt * 0.5 * (cs%PFv(i,j,k) + cs%CAv(i,j,k)))
  enddo ; enddo ; enddo
  call cpu_clock_end(id_clock_mom_update)
 
  if (cs%debug) then
    call mom_state_chksum("Predictor 2", upp, vpp, h_av, uh, vh, g, gv, us)
    call mom_accel_chksum("Predictor 2 accel", cs%CAu, cs%CAv, cs%PFu, cs%PFv, &
                          cs%diffu, cs%diffv, g, gv, us)
  endif
 
! upp <- upp + dt/2 d/dz visc d/dz upp
  call cpu_clock_begin(id_clock_vertvisc)
  call vertvisc_coef(upp, vpp, hp, forces, visc, dt*0.5, g, gv, us, cs%vertvisc_CSp, cs%OBC)
  call vertvisc(upp, vpp, hp, forces, visc, dt*0.5, cs%OBC, cs%ADp, cs%CDp, &
                g, gv, us, cs%vertvisc_CSp, waves=waves)
  call cpu_clock_end(id_clock_vertvisc)
  call pass_vector(upp, vpp, g%Domain, clock=id_clock_pass)
 
! uh = upp * hp
! h = hp + dt/2 div . uh
  call cpu_clock_begin(id_clock_continuity)
  call continuity(upp, vpp, hp, h, uh, vh, (dt*0.5), g, gv, us, cs%continuity_CSp, obc=cs%OBC)
  call cpu_clock_end(id_clock_continuity)
  call pass_var(h, g%Domain, clock=id_clock_pass)
  call pass_vector(uh, vh, g%Domain, clock=id_clock_pass)
  ! Whenever thickness changes let the diag manager know, target grids
  ! for vertical remapping may need to be regenerated.
  call diag_update_remap_grids(cs%diag)
 
  call enable_averages(0.5*dt, time_local, cs%diag)
!   Here the second half of the thickness fluxes are offered for averaging.
  if (cs%id_uh > 0) call post_data(cs%id_uh, uh, cs%diag)
  if (cs%id_vh > 0) call post_data(cs%id_vh, vh, cs%diag)
  call disable_averaging(cs%diag)
  call enable_averages(dt, time_local, cs%diag)
 
! h_av = (h + hp)/2
  do k=1,nz
    do j=js-2,je+2 ; do i=is-2,ie+2
      h_av(i,j,k) = 0.5*(h(i,j,k) + hp(i,j,k))
    enddo ; enddo
    do j=js-2,je+2 ; do i=isq-2,ieq+2
      uhtr(i,j,k) = uhtr(i,j,k) + 0.5*dt*uh(i,j,k)
    enddo ; enddo
    do j=jsq-2,jeq+2 ; do i=is-2,ie+2
      vhtr(i,j,k) = vhtr(i,j,k) + 0.5*dt*vh(i,j,k)
    enddo ; enddo
  enddo
 
! CAu = -(f+zeta(upp))/h_av vh + d/dx KE(upp)
  call cpu_clock_begin(id_clock_cor)
  call coradcalc(upp, vpp, h_av, uh, vh, cs%CAu, cs%CAv, cs%OBC, cs%ADp, &
                 g, gv, us, cs%CoriolisAdv_CSp)
  call cpu_clock_end(id_clock_cor)
 
! PFu = d/dx M(h_av,T,S)
  call cpu_clock_begin(id_clock_pres)
  call pressureforce(h_av, tv, cs%PFu, cs%PFv, g, gv, us, &
                     cs%PressureForce_CSp, cs%ALE_CSp, p_surf)
  call cpu_clock_end(id_clock_pres)
 
  if (associated(cs%OBC)) then; if (cs%OBC%update_OBC) then
    call update_obc_data(cs%OBC, g, gv, us, tv, h, cs%update_OBC_CSp, time_local)
  endif; endif
 
! u = u + dt * ( PFu + CAu )
  if (associated(cs%OBC)) then
    call open_boundary_zero_normal_flow(cs%OBC, g, cs%PFu, cs%PFv)
    call open_boundary_zero_normal_flow(cs%OBC, g, cs%CAu, cs%CAv)
  endif
  do k=1,nz ; do j=js,je ; do i=isq,ieq
    u(i,j,k) = g%mask2dCu(i,j) * (u(i,j,k) + dt * (cs%PFu(i,j,k) + cs%CAu(i,j,k)))
  enddo ; enddo ; enddo
  do k=1,nz ; do j=jsq,jeq ; do i=is,ie
    v(i,j,k) = g%mask2dCv(i,j) * (v(i,j,k) + dt * (cs%PFv(i,j,k) + cs%CAv(i,j,k)))
  enddo ; enddo ; enddo
 
! u <- u + dt d/dz visc d/dz u
  call cpu_clock_begin(id_clock_vertvisc)
  call vertvisc_coef(u, v, h_av, forces, visc, dt, g, gv, us, cs%vertvisc_CSp, cs%OBC)
  call vertvisc(u, v, h_av, forces, visc, dt, cs%OBC, cs%ADp, cs%CDp, &
                g, gv, us, cs%vertvisc_CSp, cs%taux_bot, cs%tauy_bot, waves=waves)
  call cpu_clock_end(id_clock_vertvisc)
  call pass_vector(u, v, g%Domain, clock=id_clock_pass)
 
  if (cs%debug) then
    call mom_state_chksum("Corrector", u, v, h, uh, vh, g, gv, us)
    call mom_accel_chksum("Corrector accel", cs%CAu, cs%CAv, cs%PFu, cs%PFv, &
                          cs%diffu, cs%diffv, g, gv, us)
  endif
 
  if (gv%Boussinesq) then
    do j=js,je ; do i=is,ie ; eta_av(i,j) = -gv%Z_to_H*g%bathyT(i,j) ; enddo ; enddo
  else
    do j=js,je ; do i=is,ie ; eta_av(i,j) = 0.0 ; enddo ; enddo
  endif
  do k=1,nz ; do j=js,je ; do i=is,ie
    eta_av(i,j) = eta_av(i,j) + h_av(i,j,k)
  enddo ; enddo ; enddo
 
  if (dyn_p_surf) deallocate(p_surf)
 
!   Here various terms used in to update the momentum equations are
! offered for averaging.
  if (cs%id_PFu > 0) call post_data(cs%id_PFu, cs%PFu, cs%diag)
  if (cs%id_PFv > 0) call post_data(cs%id_PFv, cs%PFv, cs%diag)
  if (cs%id_CAu > 0) call post_data(cs%id_CAu, cs%CAu, cs%diag)
  if (cs%id_CAv > 0) call post_data(cs%id_CAv, cs%CAv, cs%diag)