mom_dynamics_split_rk2 Module Reference

Detailed Description

Time step the adiabatic dynamic core of MOM using RK2 method.

This file time steps the adiabatic dynamic core by splitting between baroclinic and barotropic modes. It uses a pseudo-second order Runge-Kutta time stepping scheme for the baroclinic momentum equation and a forward-backward coupling between the baroclinic momentum and continuity equations. This split time-stepping scheme is described in detail in Hallberg (JCP, 1997). Additional issues related to exact tracer conservation and how to ensure consistency between the barotropic and layered estimates of the free surface height are described in Hallberg and Adcroft (Ocean Modelling, 2009). This was the time stepping code that is used for most GOLD applications, including GFDL's ESM2G Earth system model, and all of the examples provided with the MOM code (although several of these solutions are routinely verified by comparison with the slower unsplit schemes).

The subroutine step_MOM_dyn_split_RK2 actually does the time stepping, while register_restarts_dyn_split_RK2 sets the fields that are found in a full restart file with this scheme, and initialize_dyn_split_RK2 initializes the cpu clocks that are used in this module. For largely historical reasons, this module does not have its own control structure, but shares the same control structure with MOM.F90 and the other MOM_dynamics_... modules.

Data Types
type	mom_dyn_split_rk2_cs
	MOM_dynamics_split_RK2 module control structure. More...
integer	id_clock_cor
	CPU time clock IDs.
integer	id_clock_pres
	CPU time clock IDs.
integer	id_clock_vertvisc
	CPU time clock IDs.
integer	id_clock_horvisc
	CPU time clock IDs.
integer	id_clock_mom_update
	CPU time clock IDs.
integer	id_clock_continuity
	CPU time clock IDs.
integer	id_clock_thick_diff
	CPU time clock IDs.
integer	id_clock_btstep
	CPU time clock IDs.
integer	id_clock_btcalc
	CPU time clock IDs.
integer	id_clock_btforce
	CPU time clock IDs.
integer	id_clock_pass
	CPU time clock IDs.
integer	id_clock_pass_init
	CPU time clock IDs.
subroutine, public	step_mom_dyn_split_rk2 (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, calc_dtbt, VarMix, MEKE, thickness_diffuse_CSp, Waves)
	RK2 splitting for time stepping MOM adiabatic dynamics. More...
subroutine, public	register_restarts_dyn_split_rk2 (HI, GV, param_file, CS, restart_CS, uh, vh)
	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. More...
subroutine, public	initialize_dyn_split_rk2 (u, v, h, uh, vh, eta, Time, G, GV, US, param_file, diag, CS, restart_CS, dt, Accel_diag, Cont_diag, MIS, VarMix, MEKE, thickness_diffuse_CSp, OBC, update_OBC_CSp, ALE_CSp, setVisc_CSp, visc, dirs, ntrunc, calc_dtbt)
	This subroutine initializes all of the variables that are used by this dynamic core, including diagnostics and the cpu clocks. More...
subroutine, public	end_dyn_split_rk2 (CS)
	Close the dyn_split_RK2 module. More...

Function/Subroutine Documentation

end_dyn_split_rk2()

subroutine, public mom_dynamics_split_rk2::end_dyn_split_rk2

(

type(mom_dyn_split_rk2_cs), pointer

CS

)

Close the dyn_split_RK2 module.

Parameters

cs	module control structure

Definition at line 1234 of file MOM_dynamics_split_RK2.F90.

  type(MOM_dyn_split_RK2_CS), pointer :: CS  !< module control structure
 
  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)
 
  call dealloc_bt_cont_type(cs%BT_cont)
 
  deallocate(cs)

initialize_dyn_split_rk2()

subroutine, public mom_dynamics_split_rk2::initialize_dyn_split_rk2	(	real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(inout)	u,
		real, dimension( g %isd: g %ied, g %jsdb: g %jedb, g %ke), intent(inout)	v,
		real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(inout)	h,
		real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(inout), target	uh,
		real, dimension( g %isd: g %ied, g %jsdb: g %jedb, g %ke), intent(inout), target	vh,
		real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(inout)	eta,
		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_split_rk2_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(thickness_diffuse_cs), pointer	thickness_diffuse_CSp,
		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,
		logical, intent(out)	calc_dtbt
	)

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

Parameters

[in,out]	g	ocean grid structure
[in]	gv	ocean vertical grid structure
[in]	us	A dimensional unit scaling type
[in,out]	u	zonal velocity [m s-1]
[in,out]	v	merid velocity [m s-1]
[in,out]	h	layer thickness [H ~> m or kg m-2]
[in,out]	uh	zonal volume/mass transport [H m2 s-1 ~> m3 s-1 or kg s-1]
[in,out]	vh	merid volume/mass transport [H m2 s-1 ~> m3 s-1 or kg s-1]
[in,out]	eta	free surface height or column mass [H ~> m or kg m-2]
[in]	time	current model time
[in]	param_file	parameter file for parsing
[in,out]	diag	to control diagnostics
	cs	module control structure
	restart_cs	restart control structure
[in]	dt	time step [s]
[in,out]	accel_diag	points to momentum equation terms for budget analysis
[in,out]	cont_diag	points to terms in continuity equation
[in,out]	mis	"MOM6 internal state" used to pass diagnostic pointers
	varmix	points to spatially variable viscosities
	meke	points to mesoscale eddy kinetic energy fields
	thickness_diffuse_csp	Pointer to the control structure used for the isopycnal height diffusive transport.
	obc	points to OBC related fields
	update_obc_csp	points to OBC update related fields
	ale_csp	points to ALE control structure
	setvisc_csp	points to the set_visc control structure.
[in,out]	visc	vertical viscosities, bottom drag, and related
[in]	dirs	contains 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]	calc_dtbt	If true, recalculate the barotropic time step

Definition at line 967 of file MOM_dynamics_split_RK2.F90.

  type(ocean_grid_type),            intent(inout) :: G          !< ocean grid structure
  type(verticalGrid_type),          intent(in)    :: GV         !< ocean 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          !< zonal velocity [m s-1]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                                    intent(inout) :: v          !< merid velocity [m s-1]
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) , intent(inout) :: h !< layer thickness [H ~> m or kg m-2]
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                            target, intent(inout) :: uh !< zonal volume/mass transport [H m2 s-1 ~> m3 s-1 or kg s-1]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                            target, intent(inout) :: vh !< merid volume/mass transport [H m2 s-1 ~> m3 s-1 or kg s-1]
  real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: eta        !< free surface height or column mass [H ~> m or kg m-2]
  type(time_type),          target, intent(in)    :: Time       !< current model time
  type(param_file_type),            intent(in)    :: param_file !< parameter file for parsing
  type(diag_ctrl),          target, intent(inout) :: diag       !< to control diagnostics
  type(MOM_dyn_split_RK2_CS),       pointer       :: CS         !< module control structure
  type(MOM_restart_CS),             pointer       :: restart_CS !< restart control structure
  real,                             intent(in)    :: dt         !< time step [s]
  type(accel_diag_ptrs),    target, intent(inout) :: Accel_diag !< points to momentum equation terms for
                                                                !! budget analysis
  type(cont_diag_ptrs),     target, intent(inout) :: Cont_diag  !< points to terms in continuity equation
  type(ocean_internal_state),       intent(inout) :: MIS        !< "MOM6 internal state" used to pass
                                                                !! diagnostic pointers
  type(VarMix_CS),                  pointer       :: VarMix     !< points to spatially variable viscosities
  type(MEKE_type),                  pointer       :: MEKE       !< points to mesoscale eddy kinetic energy fields
!  type(Barotropic_CS),              pointer       :: Barotropic_CSp !< Pointer to the control structure for
!                                                                !! the barotropic module
  type(thickness_diffuse_CS),       pointer       :: thickness_diffuse_CSp !< Pointer to the control structure
                                                  !! used for the isopycnal height diffusive transport.
  type(ocean_OBC_type),             pointer       :: OBC        !< points to OBC related fields
  type(update_OBC_CS),              pointer       :: update_OBC_CSp !< points to OBC update related fields
  type(ALE_CS),                     pointer       :: ALE_CSp    !< points to ALE control structure
  type(set_visc_CS),                pointer       :: setVisc_CSp !< points to the set_visc control structure.
  type(vertvisc_type),              intent(inout) :: visc       !< vertical viscosities, bottom drag, and related
  type(directories),                intent(in)    :: dirs       !< contains 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).
  logical,                          intent(out)   :: calc_dtbt  !< If true, recalculate the barotropic time step
 
  ! local variables
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_tmp
  character(len=40) :: mdl = "MOM_dynamics_split_RK2" ! This module's name.
  character(len=48) :: thickness_units, flux_units, eta_rest_name
  real :: H_rescale  ! A rescaling factor for thicknesses from the representation in
                     ! a restart file to the internal representation in this run.
  real :: uH_rescale ! A rescaling factor for thickness transports from the representation in
                     ! a restart file to the internal representation in this run.
  real :: H_convert
  type(group_pass_type) :: pass_av_h_uvh
  logical :: 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_split_RK2 called with an unassociated control structure.")
  if (cs%module_is_initialized) then
    call mom_error(warning, "initialize_dyn_split_RK2 called with a control "// &
                            "structure that has already been initialized.")
    return
  endif
  cs%module_is_initialized = .true.
 
  cs%diag => diag
 
  call get_param(param_file, mdl, "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 "//&
                 "of a  2nd-order Runga-Kutta baroclinic time stepping "//&
                 "scheme (0.5) and a backward Euler scheme (1) that is "//&
                 "used for the Coriolis and inertial terms.  BE may be "//&
                 "from 0.5 to 1, but instability may occur near 0.5. "//&
                 "BE is also applicable if SPLIT is false and USE_RK2 "//&
                 "is true.", units="nondim", default=0.6)
  call get_param(param_file, mdl, "BEGW", cs%begw, &
                 "If SPLIT is true, BEGW is a number from 0 to 1 that "//&
                 "controls the extent to which the treatment of gravity "//&
                 "waves is forward-backward (0) or simulated backward "//&
                 "Euler (1).  0 is almost always used. "//&
                 "If SPLIT is false and USE_RK2 is true, BEGW can be "//&
                 "between 0 and 0.5 to damp gravity waves.", &
                 units="nondim", default=0.0)
 
  call get_param(param_file, mdl, "SPLIT_BOTTOM_STRESS", cs%split_bottom_stress, &
                 "If true, provide the bottom stress calculated by the "//&
                 "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 "//&
                 "adjustment due to the change in the barotropic velocity "//&
                 "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., debuggingparam=.true.)
  call get_param(param_file, mdl, "DEBUG_OBC", cs%debug_OBC, default=.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, 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_split_RK2 called with setVisc_CSp unassociated.")
  cs%set_visc_CSp => setvisc_csp
  call updatecfltruncationvalue(time, cs%vertvisc_CSp, &
                                activate=is_new_run(restart_cs) )
 
  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
 
  eta_rest_name = "sfc" ; if (.not.gv%Boussinesq) eta_rest_name = "p_bot"
  if (.not. query_initialized(cs%eta, trim(eta_rest_name), 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) = -gv%Z_to_H * 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
  elseif ((gv%m_to_H_restart /= 0.0) .and. (gv%m_to_H_restart /= gv%m_to_H)) then
    h_rescale = gv%m_to_H / gv%m_to_H_restart
    do j=js,je ; do i=is,ie ; cs%eta(i,j) = h_rescale * cs%eta(i,j) ; 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 barotropic_init(u, v, h, cs%eta, time, g, gv, us, param_file, diag, &
                       cs%barotropic_CSp, restart_cs, calc_dtbt, 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, us, cs%hor_visc_CSp, &
                              obc=cs%OBC, bt=cs%barotropic_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, us, cs%continuity_CSp, obc=cs%OBC)
    call pass_var(h_tmp, g%Domain, clock=id_clock_pass_init)
    cs%h_av(:,:,:) = 0.5*(h(:,:,:) + h_tmp(:,:,:))
  else
    if (.not. query_initialized(cs%h_av,"h2",restart_cs)) then
      cs%h_av(:,:,:) = h(:,:,:)
    elseif ((gv%m_to_H_restart /= 0.0) .and. (gv%m_to_H_restart /= gv%m_to_H)) then
      h_rescale = gv%m_to_H / gv%m_to_H_restart
      do k=1,nz ; do j=js,je ; do i=is,ie ; cs%h_av(i,j,k) = h_rescale * cs%h_av(i,j,k) ; enddo ; enddo ; enddo
    endif
    if ((gv%m_to_H_restart /= 0.0) .and. (gv%m_to_H_restart /= gv%m_to_H)) then
      uh_rescale = gv%m_to_H / gv%m_to_H_restart
      do k=1,nz ; do j=js,je ; do i=g%IscB,g%IecB ; uh(i,j,k) = uh_rescale * uh(i,j,k) ; enddo ; enddo ; enddo
      do k=1,nz ; do j=g%JscB,g%JecB ; do i=is,ie ; vh(i,j,k) = uh_rescale * vh(i,j,k) ; enddo ; enddo ; enddo
    endif
  endif
 
  call cpu_clock_begin(id_clock_pass_init)
  call create_group_pass(pass_av_h_uvh, cs%u_av, cs%v_av, g%Domain, halo=2)
  call create_group_pass(pass_av_h_uvh, cs%h_av, g%Domain, halo=2)
  call create_group_pass(pass_av_h_uvh, uh, vh, g%Domain, halo=2)
  call do_group_pass(pass_av_h_uvh, g%Domain)
  call cpu_clock_end(id_clock_pass_init)
 
  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)
  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)
 
  cs%id_CAu = register_diag_field('ocean_model', 'CAu', diag%axesCuL, time, &
      'Zonal Coriolis and Advective Acceleration', 'm s-2')
  cs%id_CAv = register_diag_field('ocean_model', 'CAv', diag%axesCvL, time, &
      'Meridional Coriolis and Advective Acceleration', 'm s-2')
  cs%id_PFu = register_diag_field('ocean_model', 'PFu', diag%axesCuL, time, &
      'Zonal Pressure Force Acceleration', 'm s-2')
  cs%id_PFv = register_diag_field('ocean_model', 'PFv', diag%axesCvL, time, &
      'Meridional Pressure Force Acceleration', 'm s-2')
 
  cs%id_uav = register_diag_field('ocean_model', 'uav', diag%axesCuL, time, &
      'Barotropic-step Averaged Zonal Velocity', 'm s-1')
  cs%id_vav = register_diag_field('ocean_model', 'vav', diag%axesCvL, time, &
      'Barotropic-step Averaged Meridional Velocity', 'm s-1')
 
  cs%id_u_BT_accel = register_diag_field('ocean_model', 'u_BT_accel', diag%axesCuL, time, &
    'Barotropic Anomaly Zonal Acceleration', 'm s-1')
  cs%id_v_BT_accel = register_diag_field('ocean_model', 'v_BT_accel', diag%axesCvL, time, &
    'Barotropic Anomaly Meridional Acceleration', 'm s-1')
 
  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_split_rk2()

subroutine, public mom_dynamics_split_rk2::register_restarts_dyn_split_rk2	(	type(hor_index_type), intent(in)	HI,
		type(verticalgrid_type), intent(in)	GV,
		type(param_file_type), intent(in)	param_file,
		type(mom_dyn_split_rk2_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
	)

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.

Parameters

[in]	hi	Horizontal index structure
[in]	gv	ocean vertical grid structure
[in]	param_file	parameter file
	cs	module control structure
	restart_cs	restart control structure
[in,out]	uh	zonal volume/mass transport [H m2 s-1 ~> m3 s-1 or kg s-1]
[in,out]	vh	merid volume/mass transport [H m2 s-1 ~> m3 s-1 or kg s-1]

Definition at line 886 of file MOM_dynamics_split_RK2.F90.

  type(hor_index_type),          intent(in)    :: HI         !< Horizontal index structure
  type(verticalGrid_type),       intent(in)    :: GV         !< ocean vertical grid structure
  type(param_file_type),         intent(in)    :: param_file !< parameter file
  type(MOM_dyn_split_RK2_CS),    pointer       :: CS         !< module control structure
  type(MOM_restart_CS),          pointer       :: restart_CS !< restart control structure
  real, dimension(SZIB_(HI),SZJ_(HI),SZK_(GV)), &
                         target, intent(inout) :: uh !< zonal volume/mass transport [H m2 s-1 ~> m3 s-1 or kg s-1]
  real, dimension(SZI_(HI),SZJB_(HI),SZK_(GV)), &
                         target, intent(inout) :: vh !< merid volume/mass transport [H m2 s-1 ~> m3 s-1 or kg s-1]
 
  type(vardesc)      :: vd
  character(len=40)  :: mdl = "MOM_dynamics_split_RK2" ! 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 specific to this module would be allocated.
  if (associated(cs)) then
    call mom_error(warning, "register_restarts_dyn_split_RK2 called with an associated "// &
                             "control structure.")
    return
  endif
  allocate(cs)
 
  alloc_(cs%diffu(isdb:iedb,jsd:jed,nz)) ; cs%diffu(:,:,:) = 0.0
  alloc_(cs%diffv(isd:ied,jsdb:jedb,nz)) ; cs%diffv(:,:,:) = 0.0
  alloc_(cs%CAu(isdb:iedb,jsd:jed,nz))   ; cs%CAu(:,:,:)   = 0.0
  alloc_(cs%CAv(isd:ied,jsdb:jedb,nz))   ; cs%CAv(:,:,:)   = 0.0
  alloc_(cs%PFu(isdb:iedb,jsd:jed,nz))   ; cs%PFu(:,:,:)   = 0.0
  alloc_(cs%PFv(isd:ied,jsdb:jedb,nz))   ; cs%PFv(:,:,:)   = 0.0
 
  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_H
 
  thickness_units = get_thickness_units(gv)
  flux_units = get_flux_units(gv)
 
  if (gv%Boussinesq) then
    vd = var_desc("sfc",thickness_units,"Free surface Height",'h','1')
  else
    vd = var_desc("p_bot",thickness_units,"Bottom Pressure",'h','1')
  endif
  call register_restart_field(cs%eta, vd, .false., restart_cs)
 
  vd = var_desc("u2","m s-1","Auxiliary Zonal velocity",'u','L')
  call register_restart_field(cs%u_av, vd, .false., restart_cs)
 
  vd = var_desc("v2","m s-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","m s-2","Zonal horizontal viscous acceleration",'u','L')
  call register_restart_field(cs%diffu, vd, .false., restart_cs)
 
  vd = var_desc("diffv","m s-2","Meridional horizontal viscous acceleration",'v','L')
  call register_restart_field(cs%diffv, vd, .false., restart_cs)
 
  call register_barotropic_restarts(hi, gv, param_file, cs%barotropic_CSp, &
                                    restart_cs)
 

step_mom_dyn_split_rk2()

subroutine, public mom_dynamics_split_rk2::step_mom_dyn_split_rk2	(	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(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), 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(unit_scale_type), intent(in)	US,
		type(mom_dyn_split_rk2_cs), pointer	CS,
		logical, intent(in)	calc_dtbt,
		type(varmix_cs), pointer	VarMix,
		type(meke_type), pointer	MEKE,
		type(thickness_diffuse_cs), pointer	thickness_diffuse_CSp,
		type(wave_parameters_cs), optional, pointer	Waves
	)

RK2 splitting for time stepping MOM adiabatic dynamics.

Parameters

[in,out]	g	ocean grid structure
[in]	gv	ocean vertical grid structure
[in]	us	A dimensional unit scaling type
[in,out]	u	zonal velocity [m s-1]
[in,out]	v	merid velocity [m s-1]
[in,out]	h	layer thickness [H ~> m or kg m-2]
[in]	tv	thermodynamic type
[in,out]	visc	vertical visc, bottom drag, and related
[in]	time_local	model time at end of time step
[in]	dt	time step [s]
[in]	forces	A structure with the driving mechanical forces
	p_surf_begin	surf pressure at start of this dynamic time step [Pa]
	p_surf_end	surf pressure at end of this dynamic time step [Pa]
[in,out]	uh	zonal volume/mass transport
[in,out]	vh	merid volume/mass transport
[in,out]	uhtr	accumulatated zonal volume/mass transport
[in,out]	vhtr	accumulatated merid volume/mass transport
[out]	eta_av	free surface height or column mass time averaged over time step [H ~> m or kg m-2]
	cs	module control structure
[in]	calc_dtbt	if true, recalculate barotropic time step
	varmix	specify the spatially varying viscosities
	meke	related to mesoscale eddy kinetic energy param
	thickness_diffuse_csp	Pointer to a structure containing interface height diffusivities
	waves	A pointer to a structure containing fields related to the surface wave conditions

Definition at line 237 of file MOM_dynamics_split_RK2.F90.

  type(ocean_grid_type),             intent(inout) :: G            !< ocean grid structure
  type(verticalGrid_type),           intent(in)    :: GV           !< ocean vertical grid structure
  type(unit_scale_type),             intent(in)    :: US           !< A dimensional unit scaling type
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                             target, intent(inout) :: u            !< zonal velocity [m s-1]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                             target, intent(inout) :: v            !< merid velocity [m s-1]
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  &
                                     intent(inout) :: h            !< layer thickness [H ~> m or kg m-2]
  type(thermo_var_ptrs),             intent(in)    :: tv           !< thermodynamic type
  type(vertvisc_type),               intent(inout) :: visc         !< vertical visc, bottom drag, and related
  type(time_type),                   intent(in)    :: Time_local   !< model time at end of time step
  real,                              intent(in)    :: dt           !< time step [s]
  type(mech_forcing),                intent(in)    :: forces       !< A structure with the driving mechanical forces
  real, dimension(:,:),              pointer       :: p_surf_begin !< surf pressure at start of this dynamic
                                                                   !! time step [Pa]
  real, dimension(:,:),              pointer       :: p_surf_end   !< surf pressure at end   of this dynamic
                                                                   !! time step [Pa]
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                             target, intent(inout) :: uh           !< zonal volume/mass transport
                                                                   !! [H m2 s-1 ~> m3 s-1 or kg s-1]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                             target, intent(inout) :: vh           !< merid volume/mass transport
                                                                   !! [H m2 s-1 ~> m3 s-1 or kg s-1]
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                                     intent(inout) :: uhtr         !< accumulatated zonal volume/mass transport
                                                                   !! since last tracer advection [H m2 ~> m3 or kg]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                                     intent(inout) :: vhtr         !< accumulatated merid volume/mass transport
                                                                   !! since last tracer advection [H m2 ~> m3 or kg]
  real, dimension(SZI_(G),SZJ_(G)),  intent(out)   :: eta_av       !< free surface height or column mass time
                                                                   !! averaged over time step [H ~> m or kg m-2]
  type(MOM_dyn_split_RK2_CS),        pointer       :: CS           !< module control structure
  logical,                           intent(in)    :: calc_dtbt    !< if true, recalculate barotropic time step
  type(VarMix_CS),                   pointer       :: VarMix       !< specify the spatially varying viscosities
  type(MEKE_type),                   pointer       :: MEKE         !< related to mesoscale eddy kinetic energy param
  type(thickness_diffuse_CS),        pointer       :: thickness_diffuse_CSp!< Pointer to a structure containing
                                                                    !! interface height diffusivities
  type(wave_parameters_CS), optional, pointer      :: Waves        !< A pointer to a structure containing
                                                                   !! fields related to the surface wave conditions
 
  ! local variables
  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 velocity [m s-1].
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: vp   ! Predicted meridional velocity [m s-1].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G))  :: hp   ! Predicted thickness [H ~> m or kg m-2].
 
  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 [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 [H m2 s-1 ~> 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 [H m2 s-1 ~> 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,
    ! [H ~> 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 [m s-1].
 
  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 [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(), &
    eta => null()
 
  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 [m s-1].
    v_av, & ! The meridional velocity time-averaged over a time step [m s-1].
    h_av    ! The layer thickness time-averaged over a time step [H ~> 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.
  !---For group halo pass
  logical :: showCallTree, sym
 
  integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
  integer :: cont_stencil
  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
  idt = 1.0 / dt
 
  sym=.false.;if (g%Domain%symmetric) sym=.true.  ! switch to include symmetric domain in checksums
 
  showcalltree = calltree_showquery()
  if (showcalltree) call calltree_enter("step_MOM_dyn_split_RK2(), MOM_dynamics_split_RK2.F90")
 
  !$OMP parallel do default(shared)
  do k = 1, nz
    do j=g%jsd,g%jed   ; do i=g%isdB,g%iedB ;  up(i,j,k) = 0.0 ; enddo ; enddo
    do j=g%jsdB,g%jedB ; do i=g%isd,g%ied   ;  vp(i,j,k) = 0.0 ; enddo ; enddo
    do j=g%jsd,g%jed   ; do i=g%isd,g%ied   ;  hp(i,j,k) = h(i,j,k) ; enddo ; enddo
  enddo
 
  ! Update CFL truncation value as function of time
  call updatecfltruncationvalue(time_local, cs%vertvisc_CSp)
 
  if (cs%debug) then
    call mom_state_chksum("Start predictor ", u, v, h, uh, vh, g, gv, symmetric=sym)
    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 => forces%p_surf
  endif
 
  if (associated(cs%OBC)) then
    if (cs%debug_OBC) call open_boundary_test_extern_h(g, cs%OBC, h)
 
    do k=1,nz ; do j=g%jsd,g%jed ; do i=g%IsdB,g%IedB
      u_old_rad_obc(i,j,k) = u_av(i,j,k)
    enddo ; enddo ; enddo
    do k=1,nz ; do j=g%JsdB,g%JedB ; do i=g%isd,g%ied
      v_old_rad_obc(i,j,k) = v_av(i,j,k)
    enddo ; enddo ; enddo
  endif
 
  bt_cont_bt_thick = .false.
  if (associated(cs%BT_cont)) bt_cont_bt_thick = &
    (allocated(cs%BT_cont%h_u) .and. allocated(cs%BT_cont%h_v))
 
  if (cs%split_bottom_stress) then
    taux_bot => cs%taux_bot ; tauy_bot => cs%tauy_bot
  endif
 
  !--- begin set up for group halo pass
 
  cont_stencil = continuity_stencil(cs%continuity_CSp)
  !### Apart from circle_OBCs halo for eta could be 1, but halo>=3 is required
  !### to match circle_OBCs solutions. Why?
  call cpu_clock_begin(id_clock_pass)
  call create_group_pass(cs%pass_eta, eta, g%Domain) !### , halo=1)
  call create_group_pass(cs%pass_visc_rem, cs%visc_rem_u, cs%visc_rem_v, g%Domain, &
                         to_all+scalar_pair, cgrid_ne, halo=max(1,cont_stencil))
  call create_group_pass(cs%pass_uvp, up, vp, g%Domain, halo=max(1,cont_stencil))
  call create_group_pass(cs%pass_hp_uv, hp, g%Domain, halo=2)
  call create_group_pass(cs%pass_hp_uv, u_av, v_av, g%Domain, halo=2)
  call create_group_pass(cs%pass_hp_uv, uh(:,:,:), vh(:,:,:), g%Domain, halo=2)
 
  call create_group_pass(cs%pass_uv, u, v, g%Domain, halo=max(2,cont_stencil))
  call create_group_pass(cs%pass_h, h, g%Domain, halo=max(2,cont_stencil))
  call create_group_pass(cs%pass_av_uvh, u_av, v_av, g%Domain, halo=2)
  call create_group_pass(cs%pass_av_uvh, uh(:,:,:), vh(:,:,:), g%Domain, halo=2)
  call cpu_clock_end(id_clock_pass)
  !--- end set up for group halo pass
 
 
! 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, us, cs%PressureForce_CSp, &
                     cs%ALE_CSp, p_surf, cs%pbce, cs%eta_PF)
  if (dyn_p_surf) then
    pa_to_eta = 1.0 / gv%H_to_Pa
    !$OMP parallel do default(shared)
    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 (showcalltree) call calltree_waypoint("done with PressureForce (step_MOM_dyn_split_RK2)")
 
  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) .and. cs%debug_OBC) &
    call open_boundary_zero_normal_flow(cs%OBC, g, cs%PFu, cs%PFv)
 
  if (g%nonblocking_updates) &
    call start_group_pass(cs%pass_eta, g%Domain, clock=id_clock_pass)
 
! 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, us, cs%CoriolisAdv_CSp)
  call cpu_clock_end(id_clock_cor)
  if (showcalltree) call calltree_waypoint("done with CorAdCalc (step_MOM_dyn_split_RK2)")
 
! u_bc_accel = CAu + PFu + diffu(u[n-1])
  call cpu_clock_begin(id_clock_btforce)
  !$OMP parallel do default(shared)
  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)) + us%s_to_T*cs%diffu(i,j,k)
    enddo ; enddo
    do j=jsq,jeq ; do i=is,ie
      v_bc_accel(i,j,k) = (cs%Cav(i,j,k) + cs%PFv(i,j,k)) + us%s_to_T*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 mom_accel_chksum("pre-btstep accel", cs%CAu, cs%CAv, cs%PFu, cs%PFv, &
                          cs%diffu, cs%diffv, g, gv, us, cs%pbce, u_bc_accel, v_bc_accel, &
                          symmetric=sym)
    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
 
  call cpu_clock_begin(id_clock_vertvisc)
  !$OMP parallel do default(shared)
  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
    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, forces, visc, dt, g, gv, us, &
                      cs%set_visc_CSp)
  call disable_averaging(cs%diag)
 
  if (cs%debug) then
    call uvchksum("before vertvisc: up", up, vp, g%HI, haloshift=0, symmetric=sym)
  endif
  call vertvisc_coef(up, vp, h, forces, visc, dt, g, gv, us, cs%vertvisc_CSp, cs%OBC)
  call vertvisc_remnant(visc, cs%visc_rem_u, cs%visc_rem_v, dt, g, gv, us, cs%vertvisc_CSp)
  call cpu_clock_end(id_clock_vertvisc)
  if (showcalltree) call calltree_waypoint("done with vertvisc_coef (step_MOM_dyn_split_RK2)")
 
 
  call cpu_clock_begin(id_clock_pass)
  if (g%nonblocking_updates) then
    call complete_group_pass(cs%pass_eta, g%Domain)
    call start_group_pass(cs%pass_visc_rem, g%Domain)
  else
    call do_group_pass(cs%pass_eta, g%Domain)
    call do_group_pass(cs%pass_visc_rem, g%Domain)
  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 btcalc(h, g, gv, cs%barotropic_CSp, obc=cs%OBC)
  call bt_mass_source(h, eta, .true., g, gv, cs%barotropic_CSp)
  call cpu_clock_end(id_clock_btcalc)
 
  if (g%nonblocking_updates) &
    call complete_group_pass(cs%pass_visc_rem, g%Domain, clock=id_clock_pass)
 
! u_accel_bt = layer accelerations due to barotropic solver
  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, us, &
                    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 btcalc(h, g, gv, cs%barotropic_CSp, cs%BT_cont%h_u, cs%BT_cont%h_v, &
                  obc=cs%OBC)
    endif
    if (showcalltree) call calltree_waypoint("done with continuity[BT_cont] (step_MOM_dyn_split_RK2)")
  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 set_dtbt(g, gv, us, cs%barotropic_CSp, eta, cs%pbce)
  if (showcalltree) call calltree_enter("btstep(), MOM_barotropic.F90")
  ! This is the predictor step call to btstep.
  call btstep(u, v, eta, dt, u_bc_accel, v_bc_accel, forces, 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, us, 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)
  if (showcalltree) call calltree_leave("btstep()")
  call cpu_clock_end(id_clock_btstep)
 
! up = u + dt_pred*( u_bc_accel + u_accel_bt )
  dt_pred = dt * cs%be
  call cpu_clock_begin(id_clock_mom_update)
 
  !$OMP parallel do default(shared)
  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
    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, symmetric=sym)
    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, &
                  symmetric=sym, 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, us, cs%pbce, cs%u_accel_bt, cs%v_accel_bt, symmetric=sym)
    call mom_state_chksum("Predictor 1 init", u_init, v_init, h, uh, vh, g, gv, haloshift=2, &
                          symmetric=sym)
    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)
  if (cs%debug) then
    call uvchksum("0 before vertvisc: [uv]p", up, vp, g%HI,haloshift=0, symmetric=sym)
  endif
  call vertvisc_coef(up, vp, h, forces, visc, dt_pred, g, gv, us, cs%vertvisc_CSp, &
                     cs%OBC)
  call vertvisc(up, vp, h, forces, visc, dt_pred, cs%OBC, cs%ADp, cs%CDp, g, &
                gv, us, cs%vertvisc_CSp, cs%taux_bot, cs%tauy_bot, waves=waves)
  if (showcalltree) call calltree_waypoint("done with vertvisc (step_MOM_dyn_split_RK2)")
  if (g%nonblocking_updates) then
    call cpu_clock_end(id_clock_vertvisc)
    call start_group_pass(cs%pass_uvp, g%Domain, clock=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, us, cs%vertvisc_CSp)
  call cpu_clock_end(id_clock_vertvisc)
 
  call do_group_pass(cs%pass_visc_rem, g%Domain, clock=id_clock_pass)
  if (g%nonblocking_updates) then
    call complete_group_pass(cs%pass_uvp, g%Domain, clock=id_clock_pass)
  else
    call do_group_pass(cs%pass_uvp, g%Domain, clock=id_clock_pass)
  endif
 
  ! 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, us, 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)
  if (showcalltree) call calltree_waypoint("done with continuity (step_MOM_dyn_split_RK2)")
 
  call do_group_pass(cs%pass_hp_uv, g%Domain, clock=id_clock_pass)
 
  if (associated(cs%OBC)) then
 
    if (cs%debug) &
      call uvchksum("Pre OBC avg [uv]", u_av, v_av, g%HI, haloshift=1, symmetric=sym)
 
    call radiation_open_bdry_conds(cs%OBC, u_av, u_old_rad_obc, v_av, v_old_rad_obc, g, dt_pred)
 
    if (cs%debug) &
      call uvchksum("Post OBC avg [uv]", u_av, v_av, g%HI, haloshift=1, symmetric=sym)
 
    ! These should be done with a pass that excludes uh & vh.
!   call do_group_pass(CS%pass_hp_uv, G%Domain, clock=id_clock_pass)
  endif
 
  if (g%nonblocking_updates) then
    call start_group_pass(cs%pass_av_uvh, g%Domain, clock=id_clock_pass)
  endif
 
  ! h_av = (h + hp)/2
  !$OMP parallel do default(shared)
  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 bt_mass_source(hp, eta_pred, .false., 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.
    !$OMP parallel do default(shared)
    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, us, cs%PressureForce_CSp, &
                       cs%ALE_CSp, p_surf, cs%pbce, cs%eta_PF)
    call cpu_clock_end(id_clock_pres)
    if (showcalltree) call calltree_waypoint("done with PressureForce[hp=(1-b).h+b.h] (step_MOM_dyn_split_RK2)")
  endif
 
  if (g%nonblocking_updates) &
    call complete_group_pass(cs%pass_av_uvh, g%Domain, clock=id_clock_pass)
 
  if (bt_cont_bt_thick) then
    call btcalc(h, g, gv, cs%barotropic_CSp, cs%BT_cont%h_u, cs%BT_cont%h_v, &
                obc=cs%OBC)
    if (showcalltree) call calltree_waypoint("done with btcalc[BT_cont_BT_thick] (step_MOM_dyn_split_RK2)")
  endif
 
  if (cs%debug) then
    call mom_state_chksum("Predictor ", up, vp, hp, uh, vh, g, gv, symmetric=sym)
    call uvchksum("Predictor avg [uv]", u_av, v_av, g%HI, haloshift=1, symmetric=sym)
    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, us, cs%hor_visc_CSp, &
                            obc=cs%OBC, bt=cs%barotropic_CSp)
  call cpu_clock_end(id_clock_horvisc)
  if (showcalltree) call calltree_waypoint("done with horizontal_viscosity (step_MOM_dyn_split_RK2)")
 
! 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, us, cs%CoriolisAdv_CSp)
  call cpu_clock_end(id_clock_cor)
  if (showcalltree) call calltree_waypoint("done with CorAdCalc (step_MOM_dyn_split_RK2)")
 
! 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)
  !$OMP parallel do default(shared)
  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)) + us%s_to_T*cs%diffu(i,j,k)
    enddo ; enddo
    do j=jsq,jeq ; do i=is,ie
      v_bc_accel(i,j,k) = (cs%Cav(i,j,k) + cs%PFv(i,j,k)) + us%s_to_T*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 mom_accel_chksum("corr pre-btstep accel", cs%CAu, cs%CAv, cs%PFu, cs%PFv, &
                          cs%diffu, cs%diffv, g, gv, us, cs%pbce, u_bc_accel, v_bc_accel, &
                          symmetric=sym)
    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%BT_use_layer_fluxes) then
    uh_ptr => uh ; vh_ptr => vh ; u_ptr => u_av ; v_ptr => v_av
  endif
 
  if (showcalltree) call calltree_enter("btstep(), MOM_barotropic.F90")
  ! This is the corrector step call to btstep.
  call btstep(u, v, eta, dt, u_bc_accel, v_bc_accel, forces, 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, us, 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)
  do j=js,je ; do i=is,ie ; eta(i,j) = eta_pred(i,j) ; enddo ; enddo
  call cpu_clock_end(id_clock_btstep)
  if (showcalltree) call calltree_leave("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)
  !$OMP parallel do default(shared)
  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
    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, symmetric=sym)
    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, &
                  symmetric=sym, 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, us, cs%pbce, cs%u_accel_bt, cs%v_accel_bt, &
                          symmetric=sym)
  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, forces, visc, dt, g, gv, us, cs%vertvisc_CSp, cs%OBC)
  call vertvisc(u, v, h, forces, visc, dt, cs%OBC, cs%ADp, cs%CDp, g, gv, us, &
                cs%vertvisc_CSp, cs%taux_bot, cs%tauy_bot,waves=waves)
  if (g%nonblocking_updates) then
    call cpu_clock_end(id_clock_vertvisc)
    call start_group_pass(cs%pass_uv, g%Domain, clock=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, us, cs%vertvisc_CSp)
  call cpu_clock_end(id_clock_vertvisc)
  if (showcalltree) call calltree_waypoint("done with vertvisc (step_MOM_dyn_split_RK2)")
 
! Later, h_av = (h_in + h_out)/2, but for now use h_av to store h_in.
  !$OMP parallel do default(shared)
  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 do_group_pass(cs%pass_visc_rem, g%Domain, clock=id_clock_pass)
  if (g%nonblocking_updates) then
    call complete_group_pass(cs%pass_uv, g%Domain, clock=id_clock_pass)
  else
    call do_group_pass(cs%pass_uv, g%Domain, clock=id_clock_pass)
  endif
 
  ! uh = u_av * h
  ! h  = h + dt * div . uh
  ! 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, us, &
                  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)
  call do_group_pass(cs%pass_h, 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)
  if (showcalltree) call calltree_waypoint("done with continuity (step_MOM_dyn_split_RK2)")
 
  if (g%nonblocking_updates) then
    call start_group_pass(cs%pass_av_uvh, g%Domain, clock=id_clock_pass)
  else
    call do_group_pass(cs%pass_av_uvh, g%domain, clock=id_clock_pass)
  endif
 
  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.
  !$OMP parallel do default(shared)
  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) &
    call complete_group_pass(cs%pass_av_uvh, g%Domain, clock=id_clock_pass)
 
  !$OMP parallel do default(shared)
  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
    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 time 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 time 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%debug) then
    call mom_state_chksum("Corrector ", u, v, h, uh, vh, g, gv, symmetric=sym)
    call uvchksum("Corrector avg [uv]", u_av, v_av, g%HI,haloshift=1, symmetric=sym)
    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
 
  if (showcalltree) call calltree_leave("step_MOM_dyn_split_RK2()")