Data Types
type	mom_dyn_legacy_split_cs

Functions/Subroutines
subroutine, public	step_mom_dyn_legacy_split (u, v, h, tv, visc, Time_local, dt, fluxes, p_surf_begin, p_surf_end, dt_since_flux, dt_therm, uh, vh, uhtr, vhtr, eta_av, G, GV, CS, calc_dtbt, VarMix, MEKE)

subroutine, public	adjustments_dyn_legacy_split (u, v, h, dt, G, GV, CS)

subroutine, public	register_restarts_dyn_legacy_split (HI, GV, param_file, CS, restart_CS, uh, vh)

subroutine, public	initialize_dyn_legacy_split (u, v, h, uh, vh, eta, Time, G, GV, param_file, diag, CS, restart_CS, dt, Accel_diag, Cont_diag, MIS, VarMix, MEKE, OBC, update_OBC_CSp, ALE_CSp, setVisc_CSp, visc, dirs, ntrunc)

subroutine, public	end_dyn_legacy_split (CS)

Variables
integer	id_clock_cor

integer	id_clock_pres

integer	id_clock_vertvisc

integer	id_clock_horvisc

integer	id_clock_mom_update

integer	id_clock_continuity

integer	id_clock_thick_diff

integer	id_clock_btstep

integer	id_clock_btcalc

integer	id_clock_btforce

integer	id_clock_pass

integer	id_clock_pass_init

Function/Subroutine Documentation

◆ adjustments_dyn_legacy_split()

subroutine, public mom_dynamics_legacy_split::adjustments_dyn_legacy_split	(	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,
		real, intent(in)	dt,
		type(ocean_grid_type), intent(inout)	G,
		type(verticalgrid_type), intent(in)	GV,
		type(mom_dyn_legacy_split_cs), pointer	CS
	)

Parameters

[in,out]	g	The ocean's grid structure.
[in]	gv	The ocean's vertical grid structure.
[in,out]	u	The zonal velocity, in m s-1.
[in,out]	v	The meridional velocity, in m s-1.
[in,out]	h	Layer thicknesses, in H
[in]	dt	The baroclinic dynamics time step, in s.
	cs	The control structure set up by initialize_dyn_legacy_split.

Definition at line 1044 of file MOM_dynamics_legacy_split.F90.

References id_clock_continuity.

Referenced by mom::step_mom().

   type(ocean_grid_type),          intent(inout) :: g    !< The ocean's grid structure.
   type(verticalgrid_type),        intent(in)    :: gv   !< The ocean's vertical grid structure.
   real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                                   intent(inout) :: u    !< The zonal velocity, in m s-1.
   real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                                   intent(inout) :: v    !< The meridional velocity, in m s-1.
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  &
                                   intent(inout) :: h    !< Layer thicknesses, in H
                                                         !! (usually m or kg m-2).
   real,                           intent(in)    :: dt   !< The baroclinic dynamics time step, in s.
   type(mom_dyn_legacy_split_cs),  pointer       :: cs   !< The control structure set up by
                                                         !! initialize_dyn_legacy_split.
 
 ! Arguments: u - The zonal velocity, in m s-1.
 !  (in)      v - The meridional velocity, in m s-1.
 !  (in)      h - The layer thicknesses, in m or kg m-2, depending on
 !                whether the Boussinesq approximation is made.
 !  (in)      dt - The time step in s.
 !  (in)      G - The ocean's grid structure.
 !  (in)      GV - The ocean's vertical grid structure.
 !  (in)      CS - The control structure set up by initialize_dyn_legacy_split.
 
   ! Temporary arrays to contain layer thickness fluxes in m3 s-1 or kg s-1.
   real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: uh_temp
   real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: vh_temp
   ! A temporary array to contain layer projected thicknesses in m or kg m-2.
   real, dimension(SZI_(G),SZJ_(G),SZK_(G))  :: h_temp
   integer :: i, j, k, is, ie, js, je, nz
   is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
 
   if (cs%readjust_BT_trans) then
     call cpu_clock_begin(id_clock_continuity)
     call continuity(u, v, h, h_temp, uh_temp, vh_temp, dt, g, gv, &
                     cs%continuity_CSp, obc=cs%OBC)
     call cpu_clock_end(id_clock_continuity)
 !$OMP parallel default(none) shared(is,ie,js,je,nz,CS,uh_temp,vh_temp)
 !$OMP do
     do j=js,je
       do i=is-1,ie ; cs%uhbt_in(i,j) = uh_temp(i,j,1) ; enddo
       do k=2,nz ; do i=is-1,ie
         cs%uhbt_in(i,j) = cs%uhbt_in(i,j) + uh_temp(i,j,k)
       enddo ; enddo
     enddo
 !$OMP do
     do j=js-1,je
       do i=is,ie ; cs%vhbt_in(i,j) = vh_temp(i,j,1) ; enddo
       do k=2,nz ; do i=is,ie
         cs%vhbt_in(i,j) = cs%vhbt_in(i,j) + vh_temp(i,j,k)
       enddo ; enddo
     enddo
 !$OMP end parallel
     cs%readjust_velocity = .true.
   endif
 

Here is the caller graph for this function:

◆ end_dyn_legacy_split()

subroutine, public mom_dynamics_legacy_split::end_dyn_legacy_split ( type(mom_dyn_legacy_split_cs), pointer CS )

Definition at line 1540 of file MOM_dynamics_legacy_split.F90.

   type(mom_dyn_legacy_split_cs), pointer :: cs
 !  (inout)   CS - The control structure set up by initialize_dyn_legacy_split.
 
   dealloc_(cs%diffu) ; dealloc_(cs%diffv)
   dealloc_(cs%CAu)   ; dealloc_(cs%CAv)
   dealloc_(cs%PFu)   ; dealloc_(cs%PFv)
 
   if (associated(cs%taux_bot)) deallocate(cs%taux_bot)
   if (associated(cs%tauy_bot)) deallocate(cs%tauy_bot)
   dealloc_(cs%uhbt) ; dealloc_(cs%vhbt)
   dealloc_(cs%u_accel_bt) ; dealloc_(cs%v_accel_bt)
   dealloc_(cs%visc_rem_u) ; dealloc_(cs%visc_rem_v)
 
   dealloc_(cs%eta) ; dealloc_(cs%eta_PF) ; dealloc_(cs%pbce)
   dealloc_(cs%h_av) ; dealloc_(cs%u_av) ; dealloc_(cs%v_av)
   dealloc_(cs%uhbt_in) ; dealloc_(cs%vhbt_in)
 
   call dealloc_bt_cont_type(cs%BT_cont)
 
   deallocate(cs)

◆ initialize_dyn_legacy_split()

subroutine, public mom_dynamics_legacy_split::initialize_dyn_legacy_split	(	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,
		real, dimension(szib_(g),szj_(g),szk_(g)), intent(inout), target	uh,
		real, dimension(szi_(g),szjb_(g),szk_(g)), intent(inout), target	vh,
		real, dimension(szi_(g),szj_(g)), intent(inout)	eta,
		type(time_type), intent(in), target	Time,
		type(ocean_grid_type), intent(inout)	G,
		type(verticalgrid_type), intent(in)	GV,
		type(param_file_type), intent(in)	param_file,
		type(diag_ctrl), intent(inout), target	diag,
		type(mom_dyn_legacy_split_cs), pointer	CS,
		type(mom_restart_cs), pointer	restart_CS,
		real, intent(in)	dt,
		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(varmix_cs), pointer	VarMix,
		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
	)

Parameters

[in,out]	g	The ocean's grid structure.
[in]	gv	The ocean's vertical grid structure.
[in,out]	u	The zonal velocity, in m s-1.
[in,out]	v	The meridional velocity, in m s-1.
[in,out]	h	Layer thicknesses, in H
[in,out]	uh	The zonal volume or mass transport,
[in,out]	vh	The meridional volume or mass transport,
[in,out]	eta	The free surface height or column mass,
[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_legacy_split.
	restart_cs	A pointer to the restart control structure.
[in]	dt	The baroclinic dynamics time step, in s.
[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.
	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.
	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 1214 of file MOM_dynamics_legacy_split.F90.

References id_clock_btcalc, id_clock_btforce, id_clock_btstep, id_clock_continuity, id_clock_cor, id_clock_horvisc, id_clock_mom_update, id_clock_pass, id_clock_pass_init, id_clock_pres, and id_clock_vertvisc.

   type(ocean_grid_type),         intent(inout) :: g    !< The ocean's grid structure.
   type(verticalgrid_type),       intent(in)    :: gv   !< The ocean's vertical grid structure.
   real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                                  intent(inout) :: u    !< The zonal velocity, in m s-1.
   real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                                  intent(inout) :: v    !< The meridional velocity, in m s-1.
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)) , &
                                  intent(inout) :: h    !< Layer thicknesses, in H
                                                        !! (usually m or kg m-2).
   real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                          target, intent(inout) :: uh   !< The zonal volume or mass transport,
                                                        !! in m3 s-1 or kg s-1.
   real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                          target, intent(inout) :: vh   !< The meridional volume or mass transport,
                                                        !! in m3 s-1 or kg s-1.
   real, dimension(SZI_(G),SZJ_(G)),          &
                                  intent(inout) :: eta  !< The free surface height or column mass,
                                                        !! in 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_legacy_split_cs), pointer       :: cs          !< The control structure set up by
                                                               !! initialize_dyn_legacy_split.
   type(mom_restart_cs),          pointer       :: restart_cs  !< A pointer to the restart control
                                                               !! structure.
   real,                          intent(in)    :: dt          !< The baroclinic dynamics time step,
                                                               !! in s.
   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(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(ocean_obc_type),          pointer       :: obc         !< If open boundary conditions are
                                                           !! used, this points to the ocean_OBC_type
                                                           !! that was set up in MOM_initialization.
   type(update_obc_cs),           pointer       :: update_obc_csp !< If open boundary condition
                                                               !! updates are used, this points to
                                                               !! the appropriate control structure.
   type(ale_cs),                  pointer       :: ale_csp     !< This points to the ALE control
                                                               !! structure.
   type(set_visc_cs),             pointer       :: setvisc_csp !< This points to the set_visc control
                                                               !! structure.
   type(vertvisc_type),           intent(inout) :: visc        !< A structure containing vertical
                                          !! viscosities, bottom drag viscosities, and related fields.
   type(directories),             intent(in)    :: dirs        !< A structure containing several
                                                               !! relevant directory paths.
   integer, target,               intent(inout) :: ntrunc      !< A target for the variable that
                         !! records the number of times the velocity is truncated (this should be 0).
 
 ! Arguments: u - The zonal velocity, in m s-1.
 !  (inout)   v - The meridional velocity, in m s-1.
 !  (inout)   h - The layer thicknesses, in m or kg m-2, depending on whether
 !                the Boussinesq approximation is made.
 !  (inout)   uh - The zonal volume or mass transport, in m3 s-1 or kg s-1.
 !  (inout)   vh - The meridional volume or mass transport, in m3 s-1 or kg s-1.
 !  (out)     eta - The free surface height or column mass, in m or kg m-2.
 !  (in)      Time - The current model time.
 !  (in)      G - The ocean's grid structure.
 !  (in)      param_file - A structure indicating the open file to parse for
 !                         model parameter values.
 !  (in)      diag - A structure that is used to regulate diagnostic output.
 !  (inout)   CS - The control structure set up by initialize_dyn_legacy_split.
 !  (in)      restart_CS - A pointer to the restart control structure.
 !  (inout)   Accel_diag - A set of pointers to the various accelerations in
 !                  the momentum equations, which can be used for later derived
 !                  diagnostics, like energy budgets.
 !  (inout)   Cont_diag - A structure with pointers to various terms in the
 !                   continuity equations.
 !  (inout)   MIS - The "MOM6 Internal State" structure, used to pass around
 !                  pointers to various arrays for diagnostic purposes.
 !  (in)      VarMix - A pointer to a structure with fields that specify the
 !                     spatially variable viscosities.
 !  (inout)   MEKE - A pointer to a structure containing fields related to
 !                   the Mesoscale Eddy Kinetic Energy.
 !  (in)      OBC - If open boundary conditions are used, this points to the
 !                  ocean_OBC_type that was set up in MOM_initialization.
 !  (in)      update_OBC_CSp - If open boundary condition updates are used,
 !                  this points to the appropriate control structure.
 !  (in)      ALE_CS - This points to the ALE control structure.
 !  (in)      setVisc_CSp - This points to the set_visc control structure.
 !  (inout)   visc - A structure containing vertical viscosities, bottom drag
 !                   viscosities, and related fields.
 !  (in)      dirs - A structure containing several relevant directory paths.
 !  (in)      ntrunc - A target for the variable that records the number of times
 !                     the velocity is truncated (this should be 0).
 
   !   This subroutine initializes all of the variables that are used by this
   ! dynamic core, including diagnostics and the cpu clocks.
 
   !   This subroutine initializes any variables that are specific to this time
   ! stepping scheme, including the cpu clocks.
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_tmp
   character(len=40)  :: mdl = "MOM_dynamics_legacy_split" ! This module's name.
   character(len=48) :: thickness_units, flux_units
   logical :: adiabatic, use_tides, debug_truncations
   integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz
   integer :: isdb, iedb, jsdb, jedb
   is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
   isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
   isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
 
   if (.not.associated(cs)) call mom_error(fatal, &
       "initialize_dyn_legacy_split called with an unassociated control structure.")
   if (cs%module_is_initialized) then
     call mom_error(warning, "initialize_dyn_legacy_split 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, "TIDES", use_tides, &
                  "If true, apply tidal momentum forcing.", default=.false.)
   call get_param(param_file, mdl, "BE", cs%be, &
                  "If SPLIT is true, BE determines the relative weighting \n"//&
                  "of a  2nd-order Runga-Kutta baroclinic time stepping \n"//&
                  "scheme (0.5) and a backward Euler scheme (1) that is \n"//&
                  "used for the Coriolis and inertial terms.  BE may be \n"//&
                  "from 0.5 to 1, but instability may occur near 0.5. \n"//&
                  "BE is also applicable if SPLIT is false and USE_RK2 \n"//&
                  "is true.", units="nondim", default=0.6)
   call get_param(param_file, mdl, "BEGW", cs%begw, &
                  "If SPLIT is true, BEGW is a number from 0 to 1 that \n"//&
                  "controls the extent to which the treatment of gravity \n"//&
                  "waves is forward-backward (0) or simulated backward \n"//&
                  "Euler (1).  0 is almost always used.\n"//&
                  "If SPLIT is false and USE_RK2 is true, BEGW can be \n"//&
                  "between 0 and 0.5 to damp gravity waves.", &
                  units="nondim", default=0.0)
 
   call get_param(param_file, mdl, "FLUX_BT_COUPLING", cs%flux_BT_coupling, &
                  "If true, use mass fluxes to ensure consistency between \n"//&
                  "the baroclinic and barotropic modes. This is only used \n"//&
                  "if SPLIT is true.", default=.false.)
   call get_param(param_file, mdl, "READJUST_BT_TRANS", cs%readjust_BT_trans, &
                  "If true, make a barotropic adjustment to the layer \n"//&
                  "velocities after the thermodynamic part of the step \n"//&
                  "to ensure that the interaction between the thermodynamics \n"//&
                  "and the continuity solver do not change the barotropic \n"//&
                  "transport.  This is only used if FLUX_BT_COUPLING and \n"//&
                  "SPLIT are true.", default=.false.)
   call get_param(param_file, mdl, "SPLIT_BOTTOM_STRESS", cs%split_bottom_stress, &
                  "If true, provide the bottom stress calculated by the \n"//&
                  "vertical viscosity to the barotropic solver.", default=.false.)
   call get_param(param_file, mdl, "BT_USE_LAYER_FLUXES", cs%BT_use_layer_fluxes, &
                  "If true, use the summed layered fluxes plus an \n"//&
                  "adjustment due to the change in the barotropic velocity \n"//&
                  "in the barotropic continuity equation.", default=.true.)
   call get_param(param_file, mdl, "DEBUG", cs%debug, &
                  "If true, write out verbose debugging data.", default=.false.)
   call get_param(param_file, mdl, "ADIABATIC", adiabatic, default=.false., do_not_log=.true.)
   if (.not.cs%flux_BT_coupling .or. adiabatic) cs%readjust_BT_trans = .false.
   call get_param(param_file, mdl, "DEBUG_TRUNCATIONS", debug_truncations, &
                  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
 
   alloc_(cs%uhbt(isdb:iedb,jsd:jed))   ; cs%uhbt(:,:) = 0.0
   alloc_(cs%vhbt(isd:ied,jsdb:jedb))   ; cs%vhbt(:,:) = 0.0
   alloc_(cs%visc_rem_u(isdb:iedb,jsd:jed,nz)) ; cs%visc_rem_u(:,:,:) = 0.0
   alloc_(cs%visc_rem_v(isd:ied,jsdb:jedb,nz)) ; cs%visc_rem_v(:,:,:) = 0.0
   alloc_(cs%eta_PF(isd:ied,jsd:jed))   ; cs%eta_PF(:,:) = 0.0
   alloc_(cs%pbce(isd:ied,jsd:jed,nz))  ; cs%pbce(:,:,:) = 0.0
 
   alloc_(cs%u_accel_bt(isdb:iedb,jsd:jed,nz)) ; cs%u_accel_bt(:,:,:) = 0.0
   alloc_(cs%v_accel_bt(isd:ied,jsdb:jedb,nz)) ; cs%v_accel_bt(:,:,:) = 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
   mis%pbce => cs%pbce
   mis%u_accel_bt => cs%u_accel_bt ; mis%v_accel_bt => cs%v_accel_bt
   mis%u_av => cs%u_av ; mis%v_av => cs%v_av
 
   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
 !  Accel_diag%pbce => CS%pbce
 !  Accel_diag%u_accel_bt => CS%u_accel_bt ; Accel_diag%v_accel_bt => CS%v_accel_bt
 !  Accel_diag%u_av => CS%u_av ; Accel_diag%v_av => CS%v_av
 
   call continuity_init(time, g, gv, param_file, diag, cs%continuity_CSp)
   call coriolisadv_init(time, g,  param_file, diag, cs%ADp, cs%CoriolisAdv_CSp)
   if (use_tides) call tidal_forcing_init(time, g, param_file, cs%tides_CSp)
   call pressureforce_init(time, g, gv, param_file, diag, cs%PressureForce_CSp, &
                           cs%tides_CSp)
   call hor_visc_init(time, g, param_file, diag, cs%hor_visc_CSp)
   call vertvisc_init(mis, time, g, gv, param_file, diag, cs%ADp, dirs, &
                      ntrunc, cs%vertvisc_CSp)
   if (.not.associated(setvisc_csp)) call mom_error(fatal, &
     "initialize_dyn_legacy_split 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
 
   if (.not. query_initialized(cs%eta,"sfc",restart_cs))  then
     ! Estimate eta based on the layer thicknesses - h.  With the Boussinesq
     ! approximation, eta is the free surface height anomaly, while without it
     ! eta is the mass of ocean per unit area.  eta always has the same
     ! dimensions as h, either m or kg m-3.
     !   CS%eta(:,:) = 0.0 already from initialization.
     if (gv%Boussinesq) then
       do j=js,je ; do i=is,ie ; cs%eta(i,j) = -g%bathyT(i,j) ; enddo ; enddo
     endif
     do k=1,nz ; do j=js,je ; do i=is,ie
        cs%eta(i,j) = cs%eta(i,j) + h(i,j,k)
     enddo ; enddo ; enddo
   endif
   ! Copy eta into an output array.
   do j=js,je ; do i=is,ie ; eta(i,j) = cs%eta(i,j) ; enddo ; enddo
 
   call legacy_barotropic_init(u, v, h, cs%eta, time, g, gv, param_file, diag, &
                        cs%barotropic_CSp, restart_cs, cs%BT_cont, cs%tides_CSp)
 
   if (.not. query_initialized(cs%diffu,"diffu",restart_cs) .or. &
       .not. query_initialized(cs%diffv,"diffv",restart_cs)) &
     call horizontal_viscosity(u, v, h, cs%diffu, cs%diffv, meke, varmix, &
                               g, gv, cs%hor_visc_CSp)
   if (.not. query_initialized(cs%u_av,"u2", restart_cs) .or. &
       .not. query_initialized(cs%u_av,"v2", restart_cs)) then
     cs%u_av(:,:,:) = u(:,:,:)
     cs%v_av(:,:,:) = v(:,:,:)
   endif
 ! This call is just here to initialize uh and vh.
   if (.not. query_initialized(uh,"uh",restart_cs) .or. &
       .not. query_initialized(vh,"vh",restart_cs)) then
     h_tmp(:,:,:) = h(:,:,:)
     call continuity(u, v, h, h_tmp, uh, vh, dt, g, gv, cs%continuity_CSp, obc=cs%OBC)
     call cpu_clock_begin(id_clock_pass_init)
     call pass_var(h_tmp, g%Domain)
     call cpu_clock_end(id_clock_pass_init)
     cs%h_av(:,:,:) = 0.5*(h(:,:,:) + h_tmp(:,:,:))
   else
     if (.not. query_initialized(cs%h_av,"h2",restart_cs)) &
       cs%h_av(:,:,:) = h(:,:,:)
   endif
 
   !   Determine whether there is a barotropic transport that is to be used
   ! to adjust the layers' velocities.
   cs%readjust_velocity = .false.
   if (cs%readjust_BT_trans) then
     if (query_initialized(cs%uhbt_in,"uhbt_in",restart_cs) .and. &
         query_initialized(cs%vhbt_in,"vhbt_in",restart_cs)) then
       cs%readjust_velocity = .true.
       call cpu_clock_begin(id_clock_pass_init)
       call pass_vector(cs%uhbt_in, cs%vhbt_in, g%Domain)
       call cpu_clock_end(id_clock_pass_init)
     endif
   endif
 
   call cpu_clock_begin(id_clock_pass_init)
   call pass_vector(cs%u_av,cs%v_av, g%Domain)
   call pass_var(cs%h_av, g%Domain)
   call pass_vector(uh, vh, g%Domain)
   call cpu_clock_end(id_clock_pass_init)
 
   flux_units = get_flux_units(gv)
   cs%id_uh = register_diag_field('ocean_model', 'uh', diag%axesCuL, time, &
       'Zonal Thickness Flux', flux_units)
   cs%id_vh = register_diag_field('ocean_model', 'vh', diag%axesCvL, time, &
       'Meridional Thickness Flux', flux_units)
   cs%id_CAu = register_diag_field('ocean_model', 'CAu', diag%axesCuL, time, &
       'Zonal Coriolis and Advective Acceleration', 'meter second-2')
   cs%id_CAv = register_diag_field('ocean_model', 'CAv', diag%axesCvL, time, &
       'Meridional Coriolis and Advective Acceleration', 'meter second-2')
   cs%id_PFu = register_diag_field('ocean_model', 'PFu', diag%axesCuL, time, &
       'Zonal Pressure Force Acceleration', 'meter second-2')
   cs%id_PFv = register_diag_field('ocean_model', 'PFv', diag%axesCvL, time, &
       'Meridional Pressure Force Acceleration', 'meter second-2')
 
   cs%id_uav = register_diag_field('ocean_model', 'uav', diag%axesCuL, time, &
       'Barotropic-step Averaged Zonal Velocity', 'meter second-1')
   cs%id_vav = register_diag_field('ocean_model', 'vav', diag%axesCvL, time, &
       'Barotropic-step Averaged Meridional Velocity', 'meter second-1')
 
 
   if (cs%flux_BT_coupling) then
     cs%id_du_adj = register_diag_field('ocean_model', 'du_adj', diag%axesCuL, time, &
         'Zonal velocity adjustments due to nonstandard terms', 'meter second-1')
     cs%id_dv_adj = register_diag_field('ocean_model', 'dv_adj', diag%axesCvL, time, &
         'Meridional velocity adjustments due to nonstandard terms', 'meter second-1')
     if (cs%id_du_adj > 0) call safe_alloc_ptr(cs%ADp%du_other,isdb,iedb,jsd,jed,nz)
     if (cs%id_dv_adj > 0) call safe_alloc_ptr(cs%ADp%dv_other,isd,ied,jsdb,jedb,nz)
   endif
   cs%id_u_BT_accel = register_diag_field('ocean_model', 'u_BT_accel', diag%axesCuL, time, &
     'Barotropic Anomaly Zonal Acceleration', 'meter second-1')
   cs%id_v_BT_accel = register_diag_field('ocean_model', 'v_BT_accel', diag%axesCvL, time, &
     'Barotropic Anomaly Meridional Acceleration', 'meter second-1')
 
   if (debug_truncations) then
     if (cs%flux_BT_coupling) then
       call safe_alloc_ptr(cs%ADp%du_other,isdb,iedb,jsd,jed,nz)
       call safe_alloc_ptr(cs%ADp%dv_other,isd,ied,jsdb,jedb,nz)
     endif
   endif
 
   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)
   id_clock_btcalc = cpu_clock_id('(Ocean barotropic mode calc)', grain=clock_module)
   id_clock_btstep = cpu_clock_id('(Ocean barotropic mode stepping)', grain=clock_module)
   id_clock_btforce = cpu_clock_id('(Ocean barotropic forcing calc)', grain=clock_module)
 
 

◆ register_restarts_dyn_legacy_split()

subroutine, public mom_dynamics_legacy_split::register_restarts_dyn_legacy_split	(	type(hor_index_type), intent(in)	HI,
		type(verticalgrid_type), intent(in)	GV,
		type(param_file_type), intent(in)	param_file,
		type(mom_dyn_legacy_split_cs), pointer	CS,
		type(mom_restart_cs), pointer	restart_CS,
		real, dimension(szib_(hi),szj_(hi),szk_(gv)), intent(inout), target	uh,
		real, dimension(szi_(hi),szjb_(hi),szk_(gv)), intent(inout), target	vh
	)

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_legacy_split.
	restart_cs	A pointer to the restart control structure.
[in,out]	uh	The zonal volume or mass transport,
[in,out]	vh	The meridional volume or mass

Definition at line 1103 of file MOM_dynamics_legacy_split.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_legacy_split_cs), pointer      :: cs         !< The control structure set up by
                                                             !! initialize_dyn_legacy_split.
   type(mom_restart_cs),          pointer      :: restart_cs !< A pointer to the restart
                                                             !! control structure.
   real, dimension(SZIB_(HI),SZJ_(HI),SZK_(GV)), &
                         target, intent(inout) :: uh         !< The zonal volume or mass transport,
                                                             !! in m3 s-1 or kg s-1.
   real, dimension(SZI_(HI),SZJB_(HI),SZK_(GV)), &
                         target, intent(inout) :: vh         !< The meridional volume or mass
                                                             !! transport, in m3 s-1 or kg s-1.
 !   This subroutine sets up any auxiliary restart variables that are specific
 ! to the unsplit time stepping scheme.  All variables registered here should
 ! have the ability to be recreated if they are not present in a restart file.
 
 ! Arguments: HI - A horizontal index type structure.
 !  (in)      GV - The ocean's vertical grid structure.
 !  (in)      param_file - A structure indicating the open file to parse for
 !                         model parameter values.
 !  (inout)   CS - The control structure set up by initialize_dyn_legacy_split.
 !  (in)      restart_CS - A pointer to the restart control structure.
 !  (inout)   uh - The zonal volume or mass transport, in m3 s-1 or kg s-1.
 !  (inout)   vh - The meridional volume or mass transport, in m3 s-1 or kg s-1.
 
   type(vardesc) :: vd
   character(len=40)  :: mdl = "MOM_dynamics_legacy_split" ! This module's name.
   character(len=48) :: thickness_units, flux_units
   logical :: adiabatic, flux_bt_coupling, readjust_bt_trans
   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_legacy_split called with an associated "// &
                              "control structure.")
     return
   endif
   allocate(cs)
 
   call get_param(param_file, "MOM_legacy_split", "FLUX_BT_COUPLING", flux_bt_coupling, &
        default=.false., do_not_log=.true.)
   call get_param(param_file, "MOM_legacy_split", "READJUST_BT_TRANS", readjust_bt_trans, &
        default=.false., do_not_log=.true.)
   call get_param(param_file, "MOM_legacy_split", "ADIABATIC", adiabatic, &
        default=.false., do_not_log=.true.)
   if (.not.flux_bt_coupling .or. adiabatic) readjust_bt_trans = .false.
 
   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
 
   alloc_(cs%eta(isd:ied,jsd:jed))       ; cs%eta(:,:) = 0.0
   alloc_(cs%u_av(isdb:iedb,jsd:jed,nz)) ; cs%u_av(:,:,:) = 0.0
   alloc_(cs%v_av(isd:ied,jsdb:jedb,nz)) ; cs%v_av(:,:,:) = 0.0
   alloc_(cs%h_av(isd:ied,jsd:jed,nz))   ; cs%h_av(:,:,:) = gv%Angstrom
   alloc_(cs%uhbt_in(isdb:iedb,jsd:jed)) ; cs%uhbt_in(:,:) = 0.0
   alloc_(cs%vhbt_in(isd:ied,jsdb:jedb)) ; cs%vhbt_in(:,:) = 0.0
 
   thickness_units = get_thickness_units(gv)
   flux_units = get_flux_units(gv)
 
   vd = var_desc("sfc",thickness_units,"Free surface Height",'h','1')
   call register_restart_field(cs%eta, vd, .false., restart_cs)
 
   vd = var_desc("u2","meter second-1","Auxiliary Zonal velocity",'u','L')
   call register_restart_field(cs%u_av, vd, .false., restart_cs)
 
   vd = var_desc("v2","meter second-1","Auxiliary Meridional velocity",'v','L')
   call register_restart_field(cs%v_av, vd, .false., restart_cs)
 
   vd = var_desc("h2",thickness_units,"Auxiliary Layer Thickness",'h','L')
   call register_restart_field(cs%h_av, vd, .false., restart_cs)
 
   vd = var_desc("uh",flux_units,"Zonal thickness flux",'u','L')
   call register_restart_field(uh, vd, .false., restart_cs)
 
   vd = var_desc("vh",flux_units,"Meridional thickness flux",'v','L')
   call register_restart_field(vh, vd, .false., restart_cs)
 
   vd = var_desc("diffu","meter second-2","Zonal horizontal viscous acceleration",'u','L')
   call register_restart_field(cs%diffu, vd, .false., restart_cs)
 
   vd = var_desc("diffv","meter second-2","Meridional horizontal viscous acceleration",'v','L')
   call register_restart_field(cs%diffv, vd, .false., restart_cs)
 
   call register_legacy_barotropic_restarts(hi, gv, param_file, &
            cs%barotropic_CSp, restart_cs)
 
   if (readjust_bt_trans) then
     vd = var_desc("uhbt_in",flux_units,"Final instantaneous barotropic zonal thickness flux",&
                  'u','1')
     call register_restart_field(cs%uhbt_in, vd, .false., restart_cs)
 
     vd = var_desc("vhbt_in",flux_units,"Final instantaneous barotropic meridional thickness flux",&
                  'v','1')
     call register_restart_field(cs%vhbt_in, vd, .false., restart_cs)
   endif
 

◆ step_mom_dyn_legacy_split()

subroutine, public mom_dynamics_legacy_split::step_mom_dyn_legacy_split	(	real, dimension(szib_(g),szj_(g),szk_(g)), intent(inout), target	u,
		real, dimension(szi_(g),szjb_(g),szk_(g)), intent(inout), target	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(forcing), intent(in)	fluxes,
		real, dimension(:,:), pointer	p_surf_begin,
		real, dimension(:,:), pointer	p_surf_end,
		real, intent(in)	dt_since_flux,
		real, intent(in)	dt_therm,
		real, dimension(szib_(g),szj_(g),szk_(g)), intent(inout), target	uh,
		real, dimension(szi_(g),szjb_(g),szk_(g)), intent(inout), target	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(mom_dyn_legacy_split_cs), pointer	CS,
		logical, intent(in)	calc_dtbt,
		type(varmix_cs), pointer	VarMix,
		type(meke_type), pointer	MEKE
	)

Parameters

[in,out]	g	The ocean's grid structure.
[in]	gv	The ocean's vertical grid structure.
[in,out]	u	The zonal velocity, in m s-1.
[in,out]	v	The meridional velocity, in m s-1.
[in,out]	h	Layer thicknesses, in H (usually 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 baroclinic dynamics time step, in s
[in]	fluxes	A structure containing pointers to any possible forcing fields. Unused fields have NULL ptrs.
	p_surf_begin	A pointer (perhaps NULL) to the surface pressure at the beginning of this dynamic step, in Pa.
	p_surf_end	A pointer (perhaps NULL) to the surface pressure at the end of this dynamic step, in Pa.
[in]	dt_since_flux	The elapsed time since fluxes were applied, in s.
[in]	dt_therm	The thermodynamic time step, in s.
[in,out]	uh	The zonal volume or mass transport,
[in,out]	vh	The meridional volume or mass transport,
[in,out]	uhtr	The accumulated zonal volume or mass
[in,out]	vhtr	The accumulated meridional volume or mass
[out]	eta_av	The free surface height or column mass
	cs	The control structure set up by
[in]	calc_dtbt	If true, recalculate the barotropic time step.
	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 275 of file MOM_dynamics_legacy_split.F90.

References mom_continuity::continuity(), mom_coriolisadv::coradcalc(), id_clock_btcalc, id_clock_btforce, id_clock_btstep, id_clock_continuity, id_clock_cor, id_clock_horvisc, id_clock_mom_update, id_clock_pass, id_clock_pres, id_clock_vertvisc, mom_legacy_barotropic::legacy_set_dtbt(), mom_checksum_packages::mom_accel_chksum(), mom_open_boundary::open_boundary_zero_normal_flow(), mom_pressureforce::pressureforce(), mom_set_visc::set_viscous_ml(), mom_domains::to_all, and mom_boundary_update::update_obc_data().

   type(ocean_grid_type),      intent(inout) :: g    !< The ocean's grid structure.
   type(verticalgrid_type),    intent(in)    :: gv   !< The ocean's vertical grid structure.
   real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                       target, intent(inout) :: u    !< The zonal velocity, in m s-1.
   real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                       target, intent(inout) :: v    !< The meridional velocity, in m s-1.
   real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  &
                               intent(inout) :: h   !< Layer thicknesses, in H (usually m or kg m-2).
   type(thermo_var_ptrs),      intent(in)    :: tv   !< A structure pointing to various.
                                                     !! thermodynamic variables.
   type(vertvisc_type),        intent(inout) :: visc !< A structure containing vertical viscosities,
                                                     !! bottom drag viscosities, and related fields.
   type(time_type),            intent(in)    :: time_local !< The model time at the end
                                                           !! of the time step.
   real,                       intent(in)    :: dt   !< The baroclinic dynamics time step, in s
   type(forcing),              intent(in)    :: fluxes !< A structure containing pointers to any
                                                       !! possible forcing fields.  Unused fields
                                                       !! have NULL ptrs.
   real, dimension(:,:),       pointer       :: p_surf_begin !< A pointer (perhaps NULL) to the
                                                             !! surface pressure at the beginning
                                                             !! of this dynamic step, in Pa.
   real, dimension(:,:),       pointer       :: p_surf_end !< A pointer (perhaps NULL) to the
                                                           !! surface pressure at the end of
                                                           !! this dynamic step, in Pa.
   real,                       intent(in)    :: dt_since_flux !< The elapsed time since fluxes
                                                              !! were applied, in s.
   real,                       intent(in)    :: dt_therm !< The thermodynamic time step, in s.
   real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                       target, intent(inout) :: uh !< The zonal volume or mass transport,
                                                   !! in m3 s-1 or kg s-1.
   real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                       target, intent(inout) :: vh !< The meridional volume or mass transport,
                                                   !! in m3 s-1 or kg s-1.
   real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                               intent(inout) :: uhtr !< The accumulated zonal volume or mass
                                                     !! transport since the last tracer advection,
                                                     !! in m3 or kg.
   real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                               intent(inout) :: vhtr !< The accumulated meridional volume or mass
                                                     !! transport since the last tracer advection,
                                                     !! in m3 or kg.
   real, dimension(SZI_(G),SZJ_(G)),          &
                               intent(out)   :: eta_av !< The free surface height or column mass
                                                       !! time-averaged over a time step,
                                                       !! in m or kg m-2.
   type(mom_dyn_legacy_split_cs),             &
                               pointer       :: cs !< The control structure set up by
                                                   !! initialize_dyn_legacy_split.
   logical,                    intent(in)    :: calc_dtbt !< If true, recalculate the
                                                          !! barotropic time step.
   type(varmix_cs),            pointer       :: varmix !<A pointer to a structure with fields that
                                                       !! specify the spatially variable viscosities.
   type(meke_type),            pointer       :: meke   !< A pointer to a structure containing fields
                                                     !! related to the Mesoscale Eddy Kinetic Energy.
 ! Arguments: u - The zonal velocity, in m s-1.
 !  (inout)   v - The meridional velocity, in m s-1.
 !  (inout)   h - The layer thicknesses, in m or kg m-2, depending on
 !                whether the Boussinesq approximation is made.
 !  (in)      tv - a structure pointing to various thermodynamic variables.
 !  (inout)   visc - A structure containing vertical viscosities, bottom drag
 !                   viscosities, and related fields.
 !  (in)      Time_local - The model time at the end of the time step.
 !  (in)      dt - The time step in s.
 !  (in)      fluxes - A structure containing pointers to any possible
 !                     forcing fields.  Unused fields have NULL ptrs.
 !  (in)      p_surf_begin - A pointer (perhaps NULL) to the surface pressure
 !                     at the beginning of this dynamic step, in Pa.
 !  (in)      p_surf_end - A pointer (perhaps NULL) to the surface pressure
 !                     at the end of this dynamic step, in Pa.
 !  (in)      dt_since_flux - The elapsed time since fluxes were applied, in s.
 !  (in)      dt_therm - The thermodynamic time step, in s.
 !  (inout)   uh - The zonal volume or mass transport, in m3 s-1 or kg s-1.
 !  (inout)   vh - The meridional volume or mass transport, in m3 s-1 or kg s-1.
 !  (inout)   uhtr - The accumulated zonal volume or mass transport since the last
 !                   tracer advection, in m3 or kg.
 !  (inout)   vhtr - The accumulated meridional volume or mass transport since the last
 !                   tracer advection, in m3 or kg.
 !  (out)     eta_av - The free surface height or column mass time-averaged
 !                     over a time step, in m or kg m-2.
 !  (in)      G - The ocean's grid structure.
 !  (in)      GV - The ocean's vertical grid structure.
 !  (in)      CS - The control structure set up by initialize_dyn_legacy_split.
 !  (in)      calc_dtbt - If true, recalculate the barotropic time step.
 !  (in)      VarMix - A pointer to a structure with fields that specify the
 !                     spatially variable viscosities.
 !  (inout)   MEKE - A pointer to a structure containing fields related to
 !                   the Mesoscale Eddy Kinetic Energy.
 
   real :: dt_pred   ! The time step for the predictor part of the baroclinic
                     ! time stepping.
 
   real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: &
     up   ! Predicted zonal velocitiy in m s-1.
   real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: &
     vp   ! Predicted meridional velocitiy in m s-1.
   real, dimension(SZI_(G),SZJ_(G),SZK_(G))  :: &
     hp   ! Predicted thickness in m or kg m-2 (H).
 
   real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: u_bc_accel
   real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: v_bc_accel
     ! u_bc_accel and v_bc_accel are the summed baroclinic accelerations of each
     ! layer calculated by the non-barotropic part of the model, both in m s-2.
   real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), target :: uh_in
   real, dimension(SZI_(G),SZJB_(G),SZK_(G)), target :: vh_in
     ! uh_in and vh_in are the zonal or meridional mass transports that would be
     ! obtained using the initial velocities, both in m3 s-1 or kg s-1.
   real, dimension(SZIB_(G),SZJ_(G)) :: uhbt_out
   real, dimension(SZI_(G),SZJB_(G)) :: vhbt_out
     ! uhbt_out and vhbt_out are the vertically summed transports from the
     ! barotropic solver based on its final velocities, both in m3 s-1 or kg s-1.
   real, dimension(SZI_(G),SZJ_(G)) :: eta_pred
     ! eta_pred is the predictor value of the free surface height or column mass,
     ! in m or kg m-2.
   real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), target :: u_adj
   real, dimension(SZI_(G),SZJB_(G),SZK_(G)), target :: v_adj
     ! u_adj and v_adj are the zonal or meridional velocities after u and v
     ! have been barotropically adjusted so the resulting transports match
     ! uhbt_out and vhbt_out, both in m s-1.
   real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: u_tmp
   real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: v_tmp
     ! u_tmp and v_tmp are temporary velocities used in diagnostic calculations.
 
   real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: u_old_rad_obc
   real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: v_old_rad_obc
     ! u_old_rad_OBC and v_old_rad_OBC are the starting velocities, which are
     ! saved for use in the Flather open boundary condition code, both in m s-1.
 
   real :: pa_to_eta ! A factor that converts pressures to the units of eta.
   real, pointer, dimension(:,:) :: &
     p_surf => null(), eta_pf_start => null(), &
     taux_bot => null(), tauy_bot => null(), &
     uhbt_in, vhbt_in, eta
   real, pointer, dimension(:,:,:) :: &
     uh_ptr => null(), u_ptr => null(),  vh_ptr => null(), v_ptr => null(), &
     u_init => null(), v_init => null(), & ! Pointers to u and v or u_adj and v_adj.
     u_av, & ! The zonal velocity time-averaged over a time step, in m s-1.
     v_av, & ! The meridional velocity time-averaged over a time step, in m s-1.
     h_av    ! The layer thickness time-averaged over a time step, in m or
             ! kg m-2.
   real :: idt
   logical :: dyn_p_surf
   logical :: bt_cont_bt_thick ! If true, use the BT_cont_type to estimate the
                               ! relative weightings of the layers in calculating
                               ! the barotropic accelerations.
   integer :: pid_bbl_h, pid_eta_pf, pid_eta, pid_visc
   integer :: pid_h, pid_u, pid_u_av, pid_uh, pid_uhbt_in
   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
   u_av => cs%u_av ; v_av => cs%v_av ; h_av => cs%h_av
   eta => cs%eta ; uhbt_in => cs%uhbt_in ; vhbt_in => cs%vhbt_in
   idt = 1.0 / dt
 
   up(:,:,:) = 0.0 ; vp(:,:,:) = 0.0 ; hp(:,:,:) = h(:,:,:)
 
   if (cs%debug) then
     call mom_state_chksum("Start predictor ", u, v, h, uh, vh, g, gv)
     call check_redundant("Start predictor u ", u, v, g)
     call check_redundant("Start predictor uh ", uh, vh, g)
   endif
 
   dyn_p_surf = associated(p_surf_begin) .and. associated(p_surf_end)
   if (dyn_p_surf) then
     p_surf => p_surf_end
     call safe_alloc_ptr(eta_pf_start,g%isd,g%ied,g%jsd,g%jed)
     eta_pf_start(:,:) = 0.0
   else
     p_surf => fluxes%p_surf
   endif
 
   if (associated(cs%OBC)) then
     do k=1,nz ; do j=js,je ; do i=is-2,ie+1
       u_old_rad_obc(i,j,k) = u(i,j,k)
     enddo ; enddo ; enddo
     do k=1,nz ; do j=js-2,je+1 ; do i=is,ie
       v_old_rad_obc(i,j,k) = v(i,j,k)
     enddo ; enddo ; enddo
   endif
 
   if (ASSOCIATED(cs%ADp%du_other)) cs%ADp%du_other(:,:,:) = 0.0
   if (ASSOCIATED(cs%ADp%dv_other)) cs%ADp%dv_other(:,:,:) = 0.0
 
   bt_cont_bt_thick = .false.
   if (associated(cs%BT_cont)) bt_cont_bt_thick = &
     (associated(cs%BT_cont%h_u) .and. associated(cs%BT_cont%h_v))
 
   if (cs%split_bottom_stress) then
     taux_bot => cs%taux_bot ; tauy_bot => cs%tauy_bot
   endif
 
 ! PFu = d/dx M(h,T,S)
 ! pbce = dM/deta
   if (cs%begw == 0.0) call enable_averaging(dt, time_local, cs%diag)
   call cpu_clock_begin(id_clock_pres)
   call pressureforce(h, tv, cs%PFu, cs%PFv, g, gv, cs%PressureForce_CSp, &
                      cs%ALE_CSp, p_surf, cs%pbce, cs%eta_PF)
   if (dyn_p_surf) then
     if (gv%Boussinesq) then
       pa_to_eta = 1.0 / (gv%Rho0*gv%g_Earth)
     else
       pa_to_eta = 1.0 / gv%H_to_Pa
     endif
     do j=jsq,jeq+1 ; do i=isq,ieq+1
       eta_pf_start(i,j) = cs%eta_PF(i,j) - pa_to_eta * &
                           (p_surf_begin(i,j) - p_surf_end(i,j))
     enddo ; enddo
   endif
   call cpu_clock_end(id_clock_pres)
   call disable_averaging(cs%diag)
 
   if (g%nonblocking_updates) then
     call cpu_clock_begin(id_clock_pass)
     pid_eta_pf = pass_var_start(cs%eta_PF, g%Domain)
     pid_eta = pass_var_start(eta, g%Domain)
     if (cs%readjust_velocity) &
       pid_uhbt_in = pass_vector_start(uhbt_in, vhbt_in, g%Domain)
     call cpu_clock_end(id_clock_pass)
   endif
 
   if (associated(cs%OBC)) then; if (cs%OBC%update_OBC) then
     call update_obc_data(cs%OBC, g, gv, tv, h, cs%update_OBC_CSp, time_local)
   endif; endif
 
 ! CAu = -(f+zeta_av)/h_av vh + d/dx KE_av
   call cpu_clock_begin(id_clock_cor)
   call coradcalc(u_av, v_av, h_av, uh, vh, cs%CAu, cs%CAv, cs%OBC, cs%ADp, g, gv, &
                  cs%CoriolisAdv_CSp)
   call cpu_clock_end(id_clock_cor)
 
 ! u_bc_accel = CAu + PFu + diffu(u[n-1])
   call cpu_clock_begin(id_clock_btforce)
   do k=1,nz ; do j=js,je ; do i=isq,ieq
     u_bc_accel(i,j,k) = (cs%Cau(i,j,k) + cs%PFu(i,j,k)) + cs%diffu(i,j,k)
   enddo ; enddo ; enddo
   do k=1,nz ; do j=jsq,jeq ; do i=is,ie
     v_bc_accel(i,j,k) = (cs%Cav(i,j,k) + cs%PFv(i,j,k)) + cs%diffv(i,j,k)
   enddo ; enddo ; enddo
   if (associated(cs%OBC)) then
     call open_boundary_zero_normal_flow(cs%OBC, g, u_bc_accel, v_bc_accel)
   endif
   call cpu_clock_end(id_clock_btforce)
 
   if (cs%debug) then
     call check_redundant("pre-btstep CS%Ca ", cs%Cau, cs%Cav, g)
     call check_redundant("pre-btstep CS%PF ", cs%PFu, cs%PFv, g)
     call check_redundant("pre-btstep CS%diff ", cs%diffu, cs%diffv, g)
     call check_redundant("pre-btstep u_bc_accel ", u_bc_accel, v_bc_accel, g)
   endif
 
   if (g%nonblocking_updates) then
     call cpu_clock_begin(id_clock_pass)
     call pass_var_complete(pid_eta_pf, cs%eta_PF, g%Domain)
     call pass_var_complete(pid_eta, eta, g%Domain)
     if (cs%readjust_velocity) &
       call pass_vector_complete(pid_uhbt_in, uhbt_in, vhbt_in, g%Domain)
     call cpu_clock_end(id_clock_pass)
   endif
 
   call cpu_clock_begin(id_clock_vertvisc)
   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 * u_bc_accel(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 * v_bc_accel(i,j,k))
   enddo ; enddo ;  enddo
   call enable_averaging(dt, time_local, cs%diag)
   call set_viscous_ml(u, v, h, tv, fluxes, visc, dt, g, gv, &
                       cs%set_visc_CSp)
   call disable_averaging(cs%diag)
 
   call vertvisc_coef(up, vp, h, fluxes, visc, dt, g, gv, cs%vertvisc_CSp, cs%OBC)
   call vertvisc_remnant(visc, cs%visc_rem_u, cs%visc_rem_v, dt, g, gv, cs%vertvisc_CSp)
   call cpu_clock_end(id_clock_vertvisc)
 
   call cpu_clock_begin(id_clock_pass)
   if (g%nonblocking_updates) then
     pid_visc = pass_vector_start(cs%visc_rem_u, cs%visc_rem_v, g%Domain, &
                                  to_all+scalar_pair, cgrid_ne)
   else
     call pass_var(cs%eta_PF, g%Domain, complete=.false.)
     call pass_var(eta, g%Domain)
     if (cs%readjust_velocity) call pass_vector(uhbt_in, vhbt_in, g%Domain)
     call pass_vector(cs%visc_rem_u, cs%visc_rem_v, g%Domain, &
                      to_all+scalar_pair, cgrid_ne)
   endif
   call cpu_clock_end(id_clock_pass)
 
   call cpu_clock_begin(id_clock_btcalc)
   ! Calculate the relative layer weights for determining barotropic quantities.
   if (.not.bt_cont_bt_thick) &
     call legacy_btcalc(h, g, gv, cs%barotropic_CSp)
   call legacy_bt_mass_source(h, eta, fluxes, .true., dt_therm, dt_since_flux, &
                       g, gv, cs%barotropic_CSp)
   call cpu_clock_end(id_clock_btcalc)
 
   if (g%nonblocking_updates) then
     call cpu_clock_begin(id_clock_pass)
     call pass_vector_complete(pid_visc, cs%visc_rem_u, cs%visc_rem_v, g%Domain, &
                      to_all+scalar_pair, cgrid_ne)
     call cpu_clock_end(id_clock_pass)
   endif
 
 ! u_accel_bt = layer accelerations due to barotropic solver
   if (cs%flux_BT_coupling) then
     call cpu_clock_begin(id_clock_continuity)
     if (cs%readjust_velocity) then
       ! Adjust the input velocites so that their transports match uhbt_out & vhbt_out.
       call continuity(u, v, h, hp, uh_in, vh_in, dt, g, gv, &
                       cs%continuity_CSp, uhbt_in, vhbt_in, cs%OBC, &
                       cs%visc_rem_u, cs%visc_rem_v, u_adj, v_adj, &
                       bt_cont=cs%BT_cont)
       u_init => u_adj ; v_init => v_adj
       if (ASSOCIATED(cs%ADp%du_other)) then ; do k=1,nz ; do j=js,je ; do i=isq,ieq
         cs%ADp%du_other(i,j,k) = u_adj(i,j,k) - u(i,j,k)
       enddo ; enddo ; enddo ; endif
       if (ASSOCIATED(cs%ADp%dv_other)) then ; do k=1,nz ; do j=jsq,jeq ; do i=is,ie
         cs%ADp%dv_other(i,j,k) = v_adj(i,j,k) - v(i,j,k)
       enddo ; enddo ; enddo ; endif
       cs%readjust_velocity = .false.
     else
       call continuity(u, v, h, hp, uh_in, vh_in, dt, g, gv, &
                       cs%continuity_CSp, obc=cs%OBC, bt_cont=cs%BT_cont)
 !###   call continuity(u, v, h, hp, uh_in, vh_in, dt, G, GV, &
 !###                   CS%continuity_CSp, OBC=CS%OBC, visc_rem_u=CS%visc_rem_u, &
 !###                      visc_rem_v=CS%visc_rem_v, BT_cont=CS%BT_cont)
       u_init => u ; v_init => v
     endif
     call cpu_clock_end(id_clock_continuity)
 
     if (bt_cont_bt_thick) then
       call cpu_clock_begin(id_clock_pass)
       call pass_vector(cs%BT_cont%h_u, cs%BT_cont%h_v, g%Domain, &
                        to_all+scalar_pair, cgrid_ne)
       call cpu_clock_end(id_clock_pass)
       call legacy_btcalc(h, g, gv, cs%barotropic_CSp, cs%BT_cont%h_u, cs%BT_cont%h_v)
     endif
     call cpu_clock_begin(id_clock_btstep)
     if (calc_dtbt) call legacy_set_dtbt(g, gv, cs%barotropic_CSp, eta, cs%pbce, cs%BT_cont)
     call legacy_btstep(.true., uh_in, vh_in, eta, dt, u_bc_accel, v_bc_accel, &
                 fluxes, cs%pbce, cs%eta_PF, uh, vh, cs%u_accel_bt, &
                 cs%v_accel_bt, eta_pred, cs%uhbt, cs%vhbt, g, gv, &
                 cs%barotropic_CSp, cs%visc_rem_u, cs%visc_rem_v, &
                 uhbt_out = uhbt_out, vhbt_out = vhbt_out, obc = cs%OBC, &
                 bt_cont = cs%BT_cont, eta_pf_start = eta_pf_start, &
                 taux_bot=taux_bot, tauy_bot=tauy_bot)
     call cpu_clock_end(id_clock_btstep)
   else
 
     if (associated(cs%BT_cont) .or. cs%BT_use_layer_fluxes) then
       call cpu_clock_begin(id_clock_continuity)
       call continuity(u, v, h, hp, uh_in, vh_in, dt, g, gv, &
                       cs%continuity_CSp, obc=cs%OBC, &
                       visc_rem_u=cs%visc_rem_u, visc_rem_v=cs%visc_rem_v, &
                       bt_cont=cs%BT_cont)
       call cpu_clock_end(id_clock_continuity)
       if (bt_cont_bt_thick) then
         call cpu_clock_begin(id_clock_pass)
         call pass_vector(cs%BT_cont%h_u, cs%BT_cont%h_v, g%Domain, &
                          to_all+scalar_pair, cgrid_ne)
         call cpu_clock_end(id_clock_pass)
         call legacy_btcalc(h, g, gv, cs%barotropic_CSp, cs%BT_cont%h_u, cs%BT_cont%h_v)
       endif
     endif
 
     if (cs%BT_use_layer_fluxes) then
       uh_ptr => uh_in; vh_ptr => vh_in; u_ptr => u; v_ptr => v
     endif
 
     u_init => u ; v_init => v
     call cpu_clock_begin(id_clock_btstep)
     if (calc_dtbt) call legacy_set_dtbt(g, gv, cs%barotropic_CSp, eta, cs%pbce)
     call legacy_btstep(.false., u, v, eta, dt, u_bc_accel, v_bc_accel, &
                 fluxes, cs%pbce, cs%eta_PF, u_av, v_av, cs%u_accel_bt, &
                 cs%v_accel_bt, eta_pred, cs%uhbt, cs%vhbt, g, gv, cs%barotropic_CSp,&
                 cs%visc_rem_u, cs%visc_rem_v, obc=cs%OBC, &
                 bt_cont = cs%BT_cont, eta_pf_start=eta_pf_start, &
                 taux_bot=taux_bot, tauy_bot=tauy_bot, &
                 uh0=uh_ptr, vh0=vh_ptr, u_uh0=u_ptr, v_vh0=v_ptr)
     call cpu_clock_end(id_clock_btstep)
   endif
 
 ! up = u + dt_pred*( u_bc_accel + u_accel_bt )
   dt_pred = dt * cs%be
   call cpu_clock_begin(id_clock_mom_update)
   do k=1,nz ; do j=jsq,jeq ; do i=is,ie
     vp(i,j,k) = g%mask2dCv(i,j) * (v_init(i,j,k) + dt_pred * &
                     (v_bc_accel(i,j,k) + cs%v_accel_bt(i,j,k)))
   enddo ; enddo ;  enddo
 
   do k=1,nz ; do j=js,je ; do i=isq,ieq
     up(i,j,k) = g%mask2dCu(i,j) * (u_init(i,j,k) + dt_pred  * &
                     (u_bc_accel(i,j,k) + cs%u_accel_bt(i,j,k)))
   enddo ; enddo ;  enddo
   call cpu_clock_end(id_clock_mom_update)
 
   if (cs%debug) then
     call uvchksum("Predictor 1 [uv]", up, vp, g%HI,haloshift=0)
     call hchksum(h, "Predictor 1 h", g%HI, haloshift=1, scale=gv%H_to_m)
     call uvchksum("Predictor 1 [uv]h", uh, vh, &
                   g%HI, haloshift=2, scale=gv%H_to_m)
 !   call MOM_state_chksum("Predictor 1", up, vp, h, uh, vh, G, GV, haloshift=1)
     call mom_accel_chksum("Predictor accel", cs%CAu, cs%CAv, cs%PFu, cs%PFv, &
              cs%diffu, cs%diffv, g, gv, cs%pbce, cs%u_accel_bt, cs%v_accel_bt)
     call mom_state_chksum("Predictor 1 init", u_init, v_init, h, uh, vh, g, gv, haloshift=2)
     call check_redundant("Predictor 1 up", up, vp, g)
     call check_redundant("Predictor 1 uh", uh, vh, g)
   endif
 
 ! up <- up + dt_pred d/dz visc d/dz up
 ! u_av  <- u_av  + dt_pred d/dz visc d/dz u_av
   call cpu_clock_begin(id_clock_vertvisc)
   call vertvisc_coef(up, vp, h, fluxes, visc, dt_pred, g, gv, cs%vertvisc_CSp, &
                      cs%OBC)
   call vertvisc(up, vp, h, fluxes, visc, dt_pred, cs%OBC, cs%ADp, cs%CDp, g, &
                 gv, cs%vertvisc_CSp, cs%taux_bot, cs%tauy_bot)
   if (g%nonblocking_updates) then
     call cpu_clock_end(id_clock_vertvisc) ; call cpu_clock_begin(id_clock_pass)
     pid_u = pass_vector_start(up, vp, g%Domain)
     call cpu_clock_end(id_clock_pass) ; call cpu_clock_begin(id_clock_vertvisc)
   endif
   call vertvisc_remnant(visc, cs%visc_rem_u, cs%visc_rem_v, dt_pred, g, gv, cs%vertvisc_CSp)
   call cpu_clock_end(id_clock_vertvisc)
 
   call cpu_clock_begin(id_clock_pass)
   call pass_vector(cs%visc_rem_u, cs%visc_rem_v, g%Domain, &
                    to_all+scalar_pair, cgrid_ne)
   if (g%nonblocking_updates) then
     call pass_vector_complete(pid_u, up, vp, g%Domain)
   else
     call pass_vector(up, vp, g%Domain)
   endif
   call cpu_clock_end(id_clock_pass)
 
 ! uh = u_av * h
 ! hp = h + dt * div . uh
   call cpu_clock_begin(id_clock_continuity)
   call continuity(up, vp, h, hp, uh, vh, dt, g, gv, cs%continuity_CSp, &
                   cs%uhbt, cs%vhbt, cs%OBC, cs%visc_rem_u, &
                   cs%visc_rem_v, u_av, v_av, bt_cont=cs%BT_cont)
   call cpu_clock_end(id_clock_continuity)
 
   call cpu_clock_begin(id_clock_pass)
   call pass_var(hp, g%Domain)
   if (g%nonblocking_updates) then
     pid_u_av = pass_vector_start(u_av, v_av, g%Domain)
     pid_uh = pass_vector_start(uh(:,:,:), vh(:,:,:), g%Domain)
   else
     call pass_vector(u_av, v_av, g%Domain, complete=.false.)
     call pass_vector(uh(:,:,:), vh(:,:,:), g%Domain)
   endif
   call cpu_clock_end(id_clock_pass)
 
   if (associated(cs%OBC)) then
     call radiation_open_bdry_conds(cs%OBC, u_av, u_old_rad_obc, v_av, v_old_rad_obc, g, dt)
   endif
 
 ! 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 ; enddo
 
 ! The correction phase of the time step starts here.
   call enable_averaging(dt, time_local, cs%diag)
 
 !   Calculate a revised estimate of the free-surface height correction to be
 ! used in the next call to btstep.  This call is at this point so that
 ! hp can be changed if CS%begw /= 0.
 ! eta_cor = ...                 (hidden inside CS%barotropic_CSp)
   call cpu_clock_begin(id_clock_btcalc)
   call legacy_bt_mass_source(hp, eta_pred, fluxes, .false., dt_therm, &
                       dt_since_flux+dt, g, gv, cs%barotropic_CSp)
   call cpu_clock_end(id_clock_btcalc)
 
   if (cs%begw /= 0.0) then
     ! hp <- (1-begw)*h_in + begw*hp
     ! Back up hp to the value it would have had after a time-step of
     ! begw*dt.  hp is not used again until recalculated by continuity.
     do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1
       hp(i,j,k) = (1.0-cs%begw)*h(i,j,k) + cs%begw*hp(i,j,k)
     enddo ; enddo ; enddo
 
 ! PFu = d/dx M(hp,T,S)
 ! pbce = dM/deta
     call cpu_clock_begin(id_clock_pres)
     call pressureforce(hp, tv, cs%PFu, cs%PFv, g, gv, &
                        cs%PressureForce_CSp, cs%ALE_CSp, &
                        p_surf, cs%pbce, cs%eta_PF)
     call cpu_clock_end(id_clock_pres)
     call cpu_clock_begin(id_clock_pass)
     call pass_var(cs%eta_PF, g%Domain)
     call cpu_clock_end(id_clock_pass)
   endif
 
   if (g%nonblocking_updates) then
     call cpu_clock_begin(id_clock_pass)
     call pass_vector_complete(pid_u_av, u_av, v_av, g%Domain)
     call pass_vector_complete(pid_uh, uh(:,:,:), vh(:,:,:), g%Domain)
     call cpu_clock_end(id_clock_pass)
   endif
 
   if (associated(cs%OBC)) then; if (cs%OBC%update_OBC) then
     call update_obc_data(cs%OBC, g, gv, tv, h, cs%update_OBC_CSp, time_local)
   endif; endif
 
   if (bt_cont_bt_thick) then
     call cpu_clock_begin(id_clock_pass)
     call pass_vector(cs%BT_cont%h_u, cs%BT_cont%h_v, g%Domain, &
                      to_all+scalar_pair, cgrid_ne)
     call cpu_clock_end(id_clock_pass)
     call legacy_btcalc(h, g, gv, cs%barotropic_CSp, cs%BT_cont%h_u, cs%BT_cont%h_v)
   endif
 
   if (cs%debug) then
     call mom_state_chksum("Predictor ", up, vp, hp, uh, vh, g, gv)
     call uvchksum("Predictor avg [uv]", u_av, v_av, g%HI, haloshift=1)
     call hchksum(h_av,"Predictor avg h",g%HI,haloshift=0, scale=gv%H_to_m)
   ! call MOM_state_chksum("Predictor avg ", u_av, v_av,  h_av,uh, vh, G, GV)
     call check_redundant("Predictor up ", up, vp, g)
     call check_redundant("Predictor uh ", uh, vh, g)
   endif
 
 ! diffu = horizontal viscosity terms (u_av)
   call cpu_clock_begin(id_clock_horvisc)
   call horizontal_viscosity(u_av, v_av, h_av, cs%diffu, cs%diffv, &
                             meke, varmix, g, gv, cs%hor_visc_CSp, obc=cs%OBC)
   call cpu_clock_end(id_clock_horvisc)
 
 ! CAu = -(f+zeta_av)/h_av vh + d/dx KE_av
   call cpu_clock_begin(id_clock_cor)
   call coradcalc(u_av, v_av, h_av, uh, vh, cs%CAu, cs%CAv, cs%OBC, cs%ADp, g, gv, &
                  cs%CoriolisAdv_CSp)
   call cpu_clock_end(id_clock_cor)
 
 ! Calculate the momentum forcing terms for the barotropic equations.
 
 ! u_bc_accel = CAu + PFu + diffu(u[n-1])
   call cpu_clock_begin(id_clock_btforce)
   do k=1,nz ; do j=js,je ; do i=isq,ieq
     u_bc_accel(i,j,k) = (cs%Cau(i,j,k) + cs%PFu(i,j,k)) + cs%diffu(i,j,k)
   enddo ; enddo ; enddo
   do k=1,nz ; do j=jsq,jeq ; do i=is,ie
     v_bc_accel(i,j,k) = (cs%Cav(i,j,k) + cs%PFv(i,j,k)) + cs%diffv(i,j,k)
   enddo ; enddo ; enddo
   if (associated(cs%OBC)) then
     call open_boundary_zero_normal_flow(cs%OBC, g, u_bc_accel, v_bc_accel)
   endif
   call cpu_clock_end(id_clock_btforce)
 
   if (cs%debug) then
     call check_redundant("corr pre-btstep CS%Ca ", cs%Cau, cs%Cav, g)
     call check_redundant("corr pre-btstep CS%PF ", cs%PFu, cs%PFv, g)
     call check_redundant("corr pre-btstep CS%diff ", cs%diffu, cs%diffv, g)
     call check_redundant("corr pre-btstep u_bc_accel ", u_bc_accel, v_bc_accel, g)
   endif
 
 ! u_accel_bt = layer accelerations due to barotropic solver
 ! pbce = dM/deta
   call cpu_clock_begin(id_clock_btstep)
   if (cs%flux_BT_coupling) then
     call legacy_btstep(.true., uh_in, vh_in, eta, dt, u_bc_accel, v_bc_accel, &
                 fluxes, cs%pbce, cs%eta_PF, uh, vh, cs%u_accel_bt, &
                 cs%v_accel_bt, eta, cs%uhbt, cs%vhbt, g, gv, &
                 cs%barotropic_CSp, cs%visc_rem_u, cs%visc_rem_v, etaav=eta_av, &
                 uhbt_out = uhbt_out, vhbt_out = vhbt_out, obc=cs%OBC, &
                 bt_cont = cs%BT_cont, eta_pf_start = eta_pf_start, &
                 taux_bot=taux_bot, tauy_bot=tauy_bot)
   else
     if (cs%BT_use_layer_fluxes) then
       uh_ptr => uh ; vh_ptr => vh ; u_ptr => u_av ; v_ptr => v_av
     endif
 
     call legacy_btstep(.false., u, v, eta, dt, u_bc_accel, v_bc_accel, &
                 fluxes, cs%pbce, cs%eta_PF, u_av, v_av, cs%u_accel_bt, &
                 cs%v_accel_bt, eta, cs%uhbt, cs%vhbt, g, gv, &
                 cs%barotropic_CSp, cs%visc_rem_u, cs%visc_rem_v, &
                 etaav=eta_av, obc=cs%OBC, &
                 bt_cont = cs%BT_cont, eta_pf_start=eta_pf_start, &
                 taux_bot=taux_bot, tauy_bot=tauy_bot, &
                 uh0=uh_ptr, vh0=vh_ptr, u_uh0=u_ptr, v_vh0=v_ptr)
   endif
   call cpu_clock_end(id_clock_btstep)
 
   if (cs%debug) then
     call check_redundant("u_accel_bt ", cs%u_accel_bt, cs%v_accel_bt, g)
   endif
 
 ! u = u + dt*( u_bc_accel + u_accel_bt )
   call cpu_clock_begin(id_clock_mom_update)
   do k=1,nz ; do j=js,je ; do i=isq,ieq
     u(i,j,k) = g%mask2dCu(i,j) * (u_init(i,j,k) + dt * &
                     (u_bc_accel(i,j,k) + cs%u_accel_bt(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_init(i,j,k) + dt * &
                     (v_bc_accel(i,j,k) + cs%v_accel_bt(i,j,k)))
   enddo ; enddo ; enddo
   call cpu_clock_end(id_clock_mom_update)
 
   if (cs%debug) then
     call uvchksum("Corrector 1 [uv]", u, v, g%HI, haloshift=0)
     call hchksum(h,"Corrector 1 h",g%HI,haloshift=2, scale=gv%H_to_m)
     call uvchksum("Corrector 1 [uv]h", &
                   uh, vh, g%HI, haloshift=2, scale=gv%H_to_m)
   ! call MOM_state_chksum("Corrector 1", u, v, h, uh, vh, G, GV, haloshift=1)
     call mom_accel_chksum("Corrector accel", cs%CAu, cs%CAv, cs%PFu, cs%PFv, &
              cs%diffu, cs%diffv, g, gv, cs%pbce, cs%u_accel_bt, cs%v_accel_bt)
   endif
 
 ! u <- u + dt d/dz visc d/dz u
 ! u_av <- u_av + dt d/dz visc d/dz u_av
   call cpu_clock_begin(id_clock_vertvisc)
   call vertvisc_coef(u, v, h, fluxes, visc, dt, g, gv, cs%vertvisc_CSp, cs%OBC)
   call vertvisc(u, v, h, fluxes, visc, dt, cs%OBC, cs%ADp, cs%CDp, g, gv, &
                 cs%vertvisc_CSp, cs%taux_bot, cs%tauy_bot)
   if (g%nonblocking_updates) then
     call cpu_clock_end(id_clock_vertvisc) ; call cpu_clock_begin(id_clock_pass)
     pid_u = pass_vector_start(u, v, g%Domain)
     call cpu_clock_end(id_clock_pass) ; call cpu_clock_begin(id_clock_vertvisc)
   endif
   call vertvisc_remnant(visc, cs%visc_rem_u, cs%visc_rem_v, dt, g, gv, cs%vertvisc_CSp)
   call cpu_clock_end(id_clock_vertvisc)
 
 ! Later, h_av = (h_in + h_out)/2, but for now use h_av to store h_in.
   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)
   enddo ; enddo ; enddo
 
   call cpu_clock_begin(id_clock_pass)
   call pass_vector(cs%visc_rem_u, cs%visc_rem_v, g%Domain, &
                    to_all+scalar_pair, cgrid_ne)
   if (g%nonblocking_updates) then
     call pass_vector_complete(pid_u, u, v, g%Domain)
   else
     call pass_vector(u, v, g%Domain)
   endif
   call cpu_clock_end(id_clock_pass)
 
 ! uh = u_av * h
 ! h = h + dt * div . uh
   if (cs%flux_BT_coupling) then
     ! u_av and v_av adjusted so their mass transports match uhbt and vhbt.
     ! Also, determine the values of u and v so that their transports
     ! that agree with uhbt_out and vhbt_out.
     if (ASSOCIATED(cs%ADp%du_other)) then ; do k=1,nz ; do j=js,je ; do i=isq,ieq
       u_tmp(i,j,k) = u(i,j,k)
     enddo ; enddo ; enddo ; endif
     if (ASSOCIATED(cs%ADp%dv_other)) then ; do k=1,nz ; do j=jsq,jeq ; do i=is,ie
       v_tmp(i,j,k) = v(i,j,k)
     enddo ; enddo ; enddo ; endif
     call cpu_clock_begin(id_clock_continuity)
     call continuity(u, v, h, h, uh, vh, dt, g, gv, &
                     cs%continuity_CSp, cs%uhbt, cs%vhbt, cs%OBC, &
                     cs%visc_rem_u, cs%visc_rem_v, u_av, v_av, &
                     uhbt_out, vhbt_out, u, v)
     call cpu_clock_end(id_clock_continuity)
     ! 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)
     if (g%nonblocking_updates) then
       call cpu_clock_begin(id_clock_pass)
       pid_h = pass_var_start(h, g%Domain)
       call cpu_clock_end(id_clock_pass)
     endif
     if (ASSOCIATED(cs%ADp%du_other)) then ; do k=1,nz ; do j=js,je ; do i=isq,ieq
       cs%ADp%du_other(i,j,k) = cs%ADp%du_other(i,j,k) + (u(i,j,k) - u_tmp(i,j,k))
     enddo ; enddo ; enddo ; endif
     if (ASSOCIATED(cs%ADp%dv_other)) then ; do k=1,nz ; do j=jsq,jeq ; do i=is,ie
       cs%ADp%dv_other(i,j,k) = cs%ADp%dv_other(i,j,k) + (v(i,j,k) - v_tmp(i,j,k))
     enddo ; enddo ; enddo ; endif
 
     call cpu_clock_begin(id_clock_vertvisc)
     call vertvisc_limit_vel(u, v, h_av, cs%ADp, cs%CDp, fluxes, visc, dt, g, gv, cs%vertvisc_CSp)
     if (g%nonblocking_updates) then
       call cpu_clock_end(id_clock_vertvisc) ; call cpu_clock_begin(id_clock_pass)
       pid_u = pass_vector_start(u, v, g%Domain)
       call cpu_clock_end(id_clock_pass) ; call cpu_clock_begin(id_clock_vertvisc)
     endif
     call vertvisc_limit_vel(u_av, v_av, h_av, cs%ADp, cs%CDp, fluxes, visc, dt, g, gv, cs%vertvisc_CSp)
     call cpu_clock_end(id_clock_vertvisc)
 
     call cpu_clock_begin(id_clock_pass)
     if (g%nonblocking_updates) then
       call pass_var_complete(pid_h, h, g%Domain)
       call pass_vector_complete(pid_u, u, v, g%Domain)
     else
       call pass_var(h, g%Domain)
       call pass_vector(u, v, g%Domain, complete=.false.)
     endif
     call cpu_clock_end(id_clock_pass)
   else
     ! u_av and v_av adjusted so their mass transports match uhbt and vhbt.
     call cpu_clock_begin(id_clock_continuity)
     call continuity(u, v, h, h, uh, vh, dt, g, gv, &
                     cs%continuity_CSp, cs%uhbt, cs%vhbt, cs%OBC, &
                     cs%visc_rem_u, cs%visc_rem_v, u_av, v_av)
     call cpu_clock_end(id_clock_continuity)
     ! 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 cpu_clock_begin(id_clock_pass)
     call pass_var(h, g%Domain)
     call cpu_clock_end(id_clock_pass)
   endif
 
   call cpu_clock_begin(id_clock_pass)
   if (g%nonblocking_updates) then
     pid_uh = pass_vector_start(uh(:,:,:), vh(:,:,:), g%Domain)
     pid_u_av = pass_vector_start(u_av, v_av, g%Domain)
   else
     call pass_vector(u_av, v_av, g%Domain, complete=.false.)
     call pass_vector(uh(:,:,:), vh(:,:,:), g%Domain)
   endif
   call cpu_clock_end(id_clock_pass)
 
   if (associated(cs%OBC)) then
     call radiation_open_bdry_conds(cs%OBC, u, u_old_rad_obc, v, v_old_rad_obc, g, dt)
   endif
 
 ! h_av = (h_in + h_out)/2 . Going in to this line, h_av = h_in.
   do k=1,nz ; do j=js-2,je+2 ; do i=is-2,ie+2
     h_av(i,j,k) = 0.5*(h_av(i,j,k) + h(i,j,k))
   enddo ; enddo ; enddo
 
   if (g%nonblocking_updates) then
     call cpu_clock_begin(id_clock_pass)
     call pass_vector_complete(pid_uh, uh(:,:,:), vh(:,:,:), g%Domain)
     call pass_vector_complete(pid_u_av, u_av, v_av, g%Domain)
     call cpu_clock_end(id_clock_pass)
   endif
 
   do k=1,nz ; do j=js-2,je+2 ; do i=isq-2,ieq+2
     uhtr(i,j,k) = uhtr(i,j,k) + uh(i,j,k)*dt
   enddo ; enddo ; enddo
   do k=1,nz ; do j=jsq-2,jeq+2 ; do i=is-2,ie+2
     vhtr(i,j,k) = vhtr(i,j,k) + vh(i,j,k)*dt
   enddo ; enddo ; enddo
 
 !   The time-averaged free surface height has already been set by the last
 !  call to btstep.
 
 !   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)
   if (cs%id_uav > 0) call post_data(cs%id_uav, u_av, cs%diag)
   if (cs%id_vav > 0) call post_data(cs%id_vav, v_av, cs%diag)
   if (cs%id_u_BT_accel > 0) call post_data(cs%id_u_BT_accel, cs%u_accel_bt, cs%diag)
   if (cs%id_v_BT_accel > 0) call post_data(cs%id_v_BT_accel, cs%v_accel_bt, cs%diag)
   if (cs%id_du_adj > 0) call post_data(cs%id_du_adj, cs%ADp%du_other, cs%diag)
   if (cs%id_dv_adj > 0) call post_data(cs%id_dv_adj, cs%ADp%dv_other, cs%diag)
   if (cs%debug) then
     call mom_state_chksum("Corrector ", u, v, h, uh, vh, g, gv)
     call uvchksum("Corrector avg [uv]", u_av, v_av, g%HI, haloshift=1)
     call hchksum(h_av,"Corrector avg h",g%HI,haloshift=1, scale=gv%H_to_m)
  !  call MOM_state_chksum("Corrector avg ", u_av, v_av, h_av, uh, vh, G, GV)
   endif
 