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
	Module control structure. More...

Functions/Subroutines
subroutine, public	step_mom_dyn_split_rk2 (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)
	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, 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)
	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...

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

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 1159 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(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(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
	)

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,out]	u	zonal velocity (m/s)
[in,out]	v	merid velocity (m/s)
[in,out]	h	layer thickness (m or kg/m2)
[in,out]	uh	zonal volume/mass transport (m3 s-1 or kg s-1)
[in,out]	vh	merid volume/mass transport (m3 s-1 or kg s-1)
[in,out]	eta	free surface height or column mass (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 (sec)
[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
	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).

Definition at line 927 of file MOM_dynamics_split_RK2.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           !< ocean grid structure
  type(verticalgrid_type),                   intent(in)    :: gv          !< ocean vertical grid structure
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: u           !< zonal velocity (m/s)
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: v           !< merid velocity (m/s)
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) , intent(inout) :: h           !< layer thickness (m or kg/m2)
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), target, intent(inout) :: uh  !< zonal volume/mass transport (m3 s-1 or kg s-1)
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), target, intent(inout) :: vh  !< merid volume/mass transport (m3 s-1 or kg s-1)
  real, dimension(SZI_(G),SZJ_(G)),          intent(inout) :: eta         !< free surface height or column mass (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 (sec)
  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(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).

  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
  type(group_pass_type) :: pass_h_tmp, 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 \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, "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, "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_split_RK2 called with setVisc_CSp unassociated.")
  cs%set_visc_CSp => setvisc_csp
  call updatecfltruncationvalue(time, cs%vertvisc_CSp, activate= &
            ((dirs%input_filename(1:1) == 'n') .and. &
             (len_trim(dirs%input_filename) == 1))   )

  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) * gv%m_to_H ; 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 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, obc=cs%OBC)
  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 create_group_pass(pass_h_tmp, h_tmp, g%Domain)
    call do_group_pass(pass_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

  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)
  cs%id_uh = register_diag_field('ocean_model', 'uh', diag%axesCuL, time, &
      'Zonal Thickness Flux', flux_units, y_cell_method='sum', v_extensive=.true.)
  cs%id_vh = register_diag_field('ocean_model', 'vh', diag%axesCvL, time, &
      'Meridional Thickness Flux', flux_units, x_cell_method='sum', v_extensive=.true.)

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

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

  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 (m3/s or kg/s)
[in,out]	vh	merid volume/mass transport (m3/s or kg/s)

Definition at line 853 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 (m3/s or kg/s)
  real, dimension(SZI_(HI),SZJB_(HI),SZK_(GV)), target, intent(inout) :: vh !< merid volume/mass transport (m3/s or kg/s)

  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

  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_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(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_split_rk2_cs), pointer	CS,
		logical, intent(in)	calc_dtbt,
		type(varmix_cs), pointer	VarMix,
		type(meke_type), pointer	MEKE
	)

RK2 splitting for time stepping MOM adiabatic dynamics.

Parameters

[in,out]	g	ocean grid structure
[in]	gv	ocean vertical grid structure
[in,out]	u	zonal velocity (m/s)
[in,out]	v	merid velocity (m/s)
[in,out]	h	layer thickness (m or kg/m2)
[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 (sec)
[in]	fluxes	forcing fields
	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]	dt_since_flux	elapsed time since fluxes were applied (sec)
[in]	dt_therm	thermodynamic time step (sec)
[in,out]	uh	zonal volume/mass transport (m3/s or kg/s)
[in,out]	vh	merid volume/mass transport (m3/s or kg/s)
[in,out]	uhtr	accumulatated zonal volume/mass transport since last tracer advection (m3 or kg)
[in,out]	vhtr	accumulatated merid volume/mass transport since last tracer advection (m3 or kg)
[out]	eta_av	free surface height or column mass time averaged over time step (m or kg/m2)
	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

Definition at line 205 of file MOM_dynamics_split_RK2.F90.

References mom_error_handler::calltree_enter(), mom_error_handler::calltree_leave(), mom_error_handler::calltree_waypoint(), mom_domains::complete_group_pass(), mom_continuity::continuity_stencil(), mom_coriolisadv::coradcalc(), mom_diag_mediator::diag_update_remap_grids(), mom_hor_visc::horizontal_viscosity(), 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_checksum_packages::mom_accel_chksum(), mom_open_boundary::open_boundary_zero_normal_flow(), mom_pressureforce::pressureforce(), mom_barotropic::set_dtbt(), mom_set_visc::set_viscous_ml(), mom_domains::start_group_pass(), mom_boundary_update::update_obc_data(), and mom_vert_friction::updatecfltruncationvalue().

  type(ocean_grid_type),                     intent(inout) :: g             !< ocean grid structure
  type(verticalgrid_type),                   intent(in)    :: gv            !< ocean vertical grid structure
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), target, intent(inout) :: u     !< zonal velocity (m/s)
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), target, intent(inout) :: v     !< merid velocity (m/s)
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(inout) :: h             !< layer thickness (m or kg/m2)
  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 (sec)
  type(forcing),                             intent(in)    :: fluxes        !< forcing fields
  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,                                      intent(in)    :: dt_since_flux !< elapsed time since fluxes were applied (sec)
  real,                                      intent(in)    :: dt_therm      !< thermodynamic time step (sec)
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), target, intent(inout) :: uh    !< zonal volume/mass transport (m3/s or kg/s)
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), target, intent(inout) :: vh    !< merid volume/mass transport (m3/s or kg/s)
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: uhtr          !< accumulatated zonal volume/mass transport since last tracer advection (m3 or kg)
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: vhtr          !< accumulatated merid volume/mass transport since last tracer advection (m3 or kg)
  real, dimension(SZI_(G),SZJ_(G)),          intent(out)   :: eta_av        !< free surface height or column mass time averaged over time step (m or kg/m2)
  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

  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 in m s-1.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: vp   ! Predicted meridional velocity 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_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(), &
    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, 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.
  !---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 => fluxes%p_surf
  endif

  if (associated(cs%OBC)) then
    do k=1,nz ; do j=js-1,je+1 ; do i=is-2,ie+1
      u_old_rad_obc(i,j,k) = u_av(i,j,k)
    enddo ; enddo ; enddo
    do k=1,nz ; do j=js-2,je+1 ; do i=is-1,ie+1
      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 = &
    (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

  !--- begin set up for group halo pass

  call cpu_clock_begin(id_clock_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 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, 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, tv, h, cs%update_OBC_CSp, time_local)
  endif; endif

  if (g%nonblocking_updates) then
    call cpu_clock_begin(id_clock_pass)
    call start_group_pass(cs%pass_eta, g%Domain)
    call cpu_clock_end(id_clock_pass)
  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)
  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)) + 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)) + 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, 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, fluxes, visc, dt, g, gv, &
                      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, 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)
  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, 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 complete_group_pass(cs%pass_visc_rem, g%Domain)
    call cpu_clock_end(id_clock_pass)
  endif

! 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, &
                    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, 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, &
              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)
  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, 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, 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 (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 cpu_clock_begin(id_clock_pass)
    call start_group_pass(cs%pass_uvp, 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 do_group_pass(cs%pass_visc_rem, g%Domain)
  if (g%nonblocking_updates) then
    call complete_group_pass(cs%pass_uvp, g%Domain)
  else
    call do_group_pass(cs%pass_uvp, 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)
  if (showcalltree) call calltree_waypoint("done with continuity (step_MOM_dyn_split_RK2)")


  if (associated(cs%OBC)) then
    call cpu_clock_begin(id_clock_pass)
    ! These should be done with a pass that excludes uh & vh.
    call do_group_pass(cs%pass_hp_uv, g%Domain)
    call cpu_clock_end(id_clock_pass)

    call radiation_open_bdry_conds(cs%OBC, u_av, u_old_rad_obc, v_av, v_old_rad_obc, g, dt_pred)
  endif

  call cpu_clock_begin(id_clock_pass)
  call do_group_pass(cs%pass_hp_uv, g%Domain)
  if (g%nonblocking_updates) then
    call start_group_pass(cs%pass_av_uvh, g%Domain)
  endif
  call cpu_clock_end(id_clock_pass)

  ! 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, 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.
    !$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, &
                       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) then
    call cpu_clock_begin(id_clock_pass)
    call complete_group_pass(cs%pass_av_uvh, g%Domain)
    call cpu_clock_end(id_clock_pass)
  endif

  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, cs%hor_visc_CSp, obc=cs%OBC)
  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, 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)) + 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)) + 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, 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, &
              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, &
              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, 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, 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)
    call start_group_pass(cs%pass_uv, 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)
  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 cpu_clock_begin(id_clock_pass)
  call do_group_pass(cs%pass_visc_rem, g%Domain)
  if (g%nonblocking_updates) then
    call complete_group_pass(cs%pass_uv, g%Domain)
  else
    call do_group_pass(cs%pass_uv, g%Domain)
  endif
  call cpu_clock_end(id_clock_pass)

  ! 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, &
                  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 cpu_clock_begin(id_clock_pass)
  call do_group_pass(cs%pass_h, g%Domain)
  call cpu_clock_end(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)")

  call cpu_clock_begin(id_clock_pass)
  if (g%nonblocking_updates) then
    call start_group_pass(cs%pass_av_uvh, g%Domain)
  else
    call do_group_pass(cs%pass_av_uvh, 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.
  !$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) then
    call cpu_clock_begin(id_clock_pass)
    call complete_group_pass(cs%pass_av_uvh, g%Domain)
    call cpu_clock_end(id_clock_pass)
  endif
  !$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()")

Here is the call graph for this function:

Variable Documentation

id_clock_btcalc

integer mom_dynamics_split_rk2::id_clock_btcalc

private

Definition at line 195 of file MOM_dynamics_split_RK2.F90.

Referenced by initialize_dyn_split_rk2(), and step_mom_dyn_split_rk2().

id_clock_btforce

integer mom_dynamics_split_rk2::id_clock_btforce

private

Definition at line 195 of file MOM_dynamics_split_RK2.F90.

Referenced by initialize_dyn_split_rk2(), and step_mom_dyn_split_rk2().

id_clock_btstep

integer mom_dynamics_split_rk2::id_clock_btstep

private

Definition at line 195 of file MOM_dynamics_split_RK2.F90.

Referenced by initialize_dyn_split_rk2(), and step_mom_dyn_split_rk2().

integer :: id_clock_btstep, id_clock_btcalc, id_clock_btforce

id_clock_continuity

integer mom_dynamics_split_rk2::id_clock_continuity

private

Definition at line 194 of file MOM_dynamics_split_RK2.F90.

Referenced by initialize_dyn_split_rk2(), and step_mom_dyn_split_rk2().

integer :: id_clock_continuity, id_clock_thick_diff

id_clock_cor

integer mom_dynamics_split_rk2::id_clock_cor

Definition at line 192 of file MOM_dynamics_split_RK2.F90.

Referenced by initialize_dyn_split_rk2(), and step_mom_dyn_split_rk2().

integer :: id_clock_cor, id_clock_pres, id_clock_vertvisc

id_clock_horvisc

integer mom_dynamics_split_rk2::id_clock_horvisc

private

Definition at line 193 of file MOM_dynamics_split_RK2.F90.

Referenced by initialize_dyn_split_rk2(), and step_mom_dyn_split_rk2().

integer :: id_clock_horvisc, id_clock_mom_update

id_clock_mom_update

integer mom_dynamics_split_rk2::id_clock_mom_update

private

Definition at line 193 of file MOM_dynamics_split_RK2.F90.

Referenced by initialize_dyn_split_rk2(), and step_mom_dyn_split_rk2().

id_clock_pass

integer mom_dynamics_split_rk2::id_clock_pass

private

Definition at line 196 of file MOM_dynamics_split_RK2.F90.

Referenced by initialize_dyn_split_rk2(), and step_mom_dyn_split_rk2().

integer :: id_clock_pass, id_clock_pass_init

id_clock_pass_init

integer mom_dynamics_split_rk2::id_clock_pass_init

private

Definition at line 196 of file MOM_dynamics_split_RK2.F90.

Referenced by initialize_dyn_split_rk2().

id_clock_pres

integer mom_dynamics_split_rk2::id_clock_pres

private

Definition at line 192 of file MOM_dynamics_split_RK2.F90.

Referenced by initialize_dyn_split_rk2(), and step_mom_dyn_split_rk2().

id_clock_thick_diff

integer mom_dynamics_split_rk2::id_clock_thick_diff

private

Definition at line 194 of file MOM_dynamics_split_RK2.F90.

id_clock_vertvisc

integer mom_dynamics_split_rk2::id_clock_vertvisc

private

Definition at line 192 of file MOM_dynamics_split_RK2.F90.

Referenced by initialize_dyn_split_rk2(), and step_mom_dyn_split_rk2().