MOM6
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. More...
 
integer id_clock_pres
 CPU time clock IDs. More...
 
integer id_clock_vertvisc
 CPU time clock IDs. More...
 
integer id_clock_horvisc
 CPU time clock IDs. More...
 
integer id_clock_mom_update
 CPU time clock IDs. More...
 
integer id_clock_continuity
 CPU time clock IDs. More...
 
integer id_clock_thick_diff
 CPU time clock IDs. More...
 
integer id_clock_btstep
 CPU time clock IDs. More...
 
integer id_clock_btcalc
 CPU time clock IDs. More...
 
integer id_clock_btforce
 CPU time clock IDs. More...
 
integer id_clock_pass
 CPU time clock IDs. More...
 
integer id_clock_pass_init
 CPU time clock IDs. More...
 
subroutine, public step_mom_dyn_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
csmodule control structure

Definition at line 1254 of file MOM_dynamics_split_RK2.F90.

1254  type(MOM_dyn_split_RK2_CS), pointer :: CS !< module control structure
1255 
1256  dealloc_(cs%diffu) ; dealloc_(cs%diffv)
1257  dealloc_(cs%CAu) ; dealloc_(cs%CAv)
1258  dealloc_(cs%PFu) ; dealloc_(cs%PFv)
1259 
1260  if (associated(cs%taux_bot)) deallocate(cs%taux_bot)
1261  if (associated(cs%tauy_bot)) deallocate(cs%tauy_bot)
1262  dealloc_(cs%uhbt) ; dealloc_(cs%vhbt)
1263  dealloc_(cs%u_accel_bt) ; dealloc_(cs%v_accel_bt)
1264  dealloc_(cs%visc_rem_u) ; dealloc_(cs%visc_rem_v)
1265 
1266  dealloc_(cs%eta) ; dealloc_(cs%eta_PF) ; dealloc_(cs%pbce)
1267  dealloc_(cs%h_av) ; dealloc_(cs%u_av) ; dealloc_(cs%v_av)
1268 
1269  call dealloc_bt_cont_type(cs%BT_cont)
1270 
1271  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]gocean grid structure
[in]gvocean vertical grid structure
[in]usA dimensional unit scaling type
[in,out]uzonal velocity [L T-1 ~> m s-1]
[in,out]vmerid velocity [L T-1 ~> m s-1]
[in,out]hlayer thickness [H ~> m or kg m-2]
[in,out]uhzonal volume/mass transport [H L2 T-1 ~> m3 s-1 or kg s-1]
[in,out]vhmerid volume/mass transport [H L2 T-1 ~> m3 s-1 or kg s-1]
[in,out]etafree surface height or column mass [H ~> m or kg m-2]
[in]timecurrent model time
[in]param_fileparameter file for parsing
[in,out]diagto control diagnostics
csmodule control structure
restart_csrestart control structure
[in]dttime step [T ~> s]
[in,out]accel_diagpoints to momentum equation terms for budget analysis
[in,out]cont_diagpoints to terms in continuity equation
[in,out]mis"MOM6 internal state" used to pass diagnostic pointers
varmixpoints to spatially variable viscosities
mekepoints to mesoscale eddy kinetic energy fields
thickness_diffuse_cspPointer to the control structure used for the isopycnal height diffusive transport.
obcpoints to OBC related fields
update_obc_csppoints to OBC update related fields
ale_csppoints to ALE control structure
setvisc_csppoints to the set_visc control structure.
[in,out]viscvertical viscosities, bottom drag, and related
[in]dirscontains directory paths
[in,out]ntruncA target for the variable that records the number of times the velocity is truncated (this should be 0).
[out]calc_dtbtIf true, recalculate the barotropic time step

Definition at line 958 of file MOM_dynamics_split_RK2.F90.

958  type(ocean_grid_type), intent(inout) :: G !< ocean grid structure
959  type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
960  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
961  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
962  intent(inout) :: u !< zonal velocity [L T-1 ~> m s-1]
963  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
964  intent(inout) :: v !< merid velocity [L T-1 ~> m s-1]
965  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) , intent(inout) :: h !< layer thickness [H ~> m or kg m-2]
966  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
967  target, intent(inout) :: uh !< zonal volume/mass transport [H L2 T-1 ~> m3 s-1 or kg s-1]
968  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
969  target, intent(inout) :: vh !< merid volume/mass transport [H L2 T-1 ~> m3 s-1 or kg s-1]
970  real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: eta !< free surface height or column mass [H ~> m or kg m-2]
971  type(time_type), target, intent(in) :: Time !< current model time
972  type(param_file_type), intent(in) :: param_file !< parameter file for parsing
973  type(diag_ctrl), target, intent(inout) :: diag !< to control diagnostics
974  type(MOM_dyn_split_RK2_CS), pointer :: CS !< module control structure
975  type(MOM_restart_CS), pointer :: restart_CS !< restart control structure
976  real, intent(in) :: dt !< time step [T ~> s]
977  type(accel_diag_ptrs), target, intent(inout) :: Accel_diag !< points to momentum equation terms for
978  !! budget analysis
979  type(cont_diag_ptrs), target, intent(inout) :: Cont_diag !< points to terms in continuity equation
980  type(ocean_internal_state), intent(inout) :: MIS !< "MOM6 internal state" used to pass
981  !! diagnostic pointers
982  type(VarMix_CS), pointer :: VarMix !< points to spatially variable viscosities
983  type(MEKE_type), pointer :: MEKE !< points to mesoscale eddy kinetic energy fields
984 ! type(Barotropic_CS), pointer :: Barotropic_CSp !< Pointer to the control structure for
985 ! !! the barotropic module
986  type(thickness_diffuse_CS), pointer :: thickness_diffuse_CSp !< Pointer to the control structure
987  !! used for the isopycnal height diffusive transport.
988  type(ocean_OBC_type), pointer :: OBC !< points to OBC related fields
989  type(update_OBC_CS), pointer :: update_OBC_CSp !< points to OBC update related fields
990  type(ALE_CS), pointer :: ALE_CSp !< points to ALE control structure
991  type(set_visc_CS), pointer :: setVisc_CSp !< points to the set_visc control structure.
992  type(vertvisc_type), intent(inout) :: visc !< vertical viscosities, bottom drag, and related
993  type(directories), intent(in) :: dirs !< contains directory paths
994  integer, target, intent(inout) :: ntrunc !< A target for the variable that records
995  !! the number of times the velocity is
996  !! truncated (this should be 0).
997  logical, intent(out) :: calc_dtbt !< If true, recalculate the barotropic time step
998 
999  ! local variables
1000  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_tmp
1001  character(len=40) :: mdl = "MOM_dynamics_split_RK2" ! This module's name.
1002  character(len=48) :: thickness_units, flux_units, eta_rest_name
1003  real :: H_rescale ! A rescaling factor for thicknesses from the representation in
1004  ! a restart file to the internal representation in this run.
1005  real :: vel_rescale ! A rescaling factor for velocities from the representation in
1006  ! a restart file to the internal representation in this run.
1007  real :: uH_rescale ! A rescaling factor for thickness transports from the representation in
1008  ! a restart file to the internal representation in this run.
1009  real :: accel_rescale ! A rescaling factor for accelerations from the representation in
1010  ! a restart file to the internal representation in this run.
1011  real :: H_convert
1012  type(group_pass_type) :: pass_av_h_uvh
1013  logical :: use_tides, debug_truncations
1014 
1015  integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz
1016  integer :: IsdB, IedB, JsdB, JedB
1017  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
1018  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
1019  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
1020 
1021  if (.not.associated(cs)) call mom_error(fatal, &
1022  "initialize_dyn_split_RK2 called with an unassociated control structure.")
1023  if (cs%module_is_initialized) then
1024  call mom_error(warning, "initialize_dyn_split_RK2 called with a control "// &
1025  "structure that has already been initialized.")
1026  return
1027  endif
1028  cs%module_is_initialized = .true.
1029 
1030  cs%diag => diag
1031 
1032  call get_param(param_file, mdl, "TIDES", use_tides, &
1033  "If true, apply tidal momentum forcing.", default=.false.)
1034  call get_param(param_file, mdl, "BE", cs%be, &
1035  "If SPLIT is true, BE determines the relative weighting "//&
1036  "of a 2nd-order Runga-Kutta baroclinic time stepping "//&
1037  "scheme (0.5) and a backward Euler scheme (1) that is "//&
1038  "used for the Coriolis and inertial terms. BE may be "//&
1039  "from 0.5 to 1, but instability may occur near 0.5. "//&
1040  "BE is also applicable if SPLIT is false and USE_RK2 "//&
1041  "is true.", units="nondim", default=0.6)
1042  call get_param(param_file, mdl, "BEGW", cs%begw, &
1043  "If SPLIT is true, BEGW is a number from 0 to 1 that "//&
1044  "controls the extent to which the treatment of gravity "//&
1045  "waves is forward-backward (0) or simulated backward "//&
1046  "Euler (1). 0 is almost always used. "//&
1047  "If SPLIT is false and USE_RK2 is true, BEGW can be "//&
1048  "between 0 and 0.5 to damp gravity waves.", &
1049  units="nondim", default=0.0)
1050 
1051  call get_param(param_file, mdl, "SPLIT_BOTTOM_STRESS", cs%split_bottom_stress, &
1052  "If true, provide the bottom stress calculated by the "//&
1053  "vertical viscosity to the barotropic solver.", default=.false.)
1054  call get_param(param_file, mdl, "BT_USE_LAYER_FLUXES", cs%BT_use_layer_fluxes, &
1055  "If true, use the summed layered fluxes plus an "//&
1056  "adjustment due to the change in the barotropic velocity "//&
1057  "in the barotropic continuity equation.", default=.true.)
1058  call get_param(param_file, mdl, "DEBUG", cs%debug, &
1059  "If true, write out verbose debugging data.", &
1060  default=.false., debuggingparam=.true.)
1061  call get_param(param_file, mdl, "DEBUG_OBC", cs%debug_OBC, default=.false.)
1062  call get_param(param_file, mdl, "DEBUG_TRUNCATIONS", debug_truncations, &
1063  default=.false.)
1064 
1065  allocate(cs%taux_bot(isdb:iedb,jsd:jed)) ; cs%taux_bot(:,:) = 0.0
1066  allocate(cs%tauy_bot(isd:ied,jsdb:jedb)) ; cs%tauy_bot(:,:) = 0.0
1067 
1068  alloc_(cs%uhbt(isdb:iedb,jsd:jed)) ; cs%uhbt(:,:) = 0.0
1069  alloc_(cs%vhbt(isd:ied,jsdb:jedb)) ; cs%vhbt(:,:) = 0.0
1070  alloc_(cs%visc_rem_u(isdb:iedb,jsd:jed,nz)) ; cs%visc_rem_u(:,:,:) = 0.0
1071  alloc_(cs%visc_rem_v(isd:ied,jsdb:jedb,nz)) ; cs%visc_rem_v(:,:,:) = 0.0
1072  alloc_(cs%eta_PF(isd:ied,jsd:jed)) ; cs%eta_PF(:,:) = 0.0
1073  alloc_(cs%pbce(isd:ied,jsd:jed,nz)) ; cs%pbce(:,:,:) = 0.0
1074 
1075  alloc_(cs%u_accel_bt(isdb:iedb,jsd:jed,nz)) ; cs%u_accel_bt(:,:,:) = 0.0
1076  alloc_(cs%v_accel_bt(isd:ied,jsdb:jedb,nz)) ; cs%v_accel_bt(:,:,:) = 0.0
1077 
1078  mis%diffu => cs%diffu
1079  mis%diffv => cs%diffv
1080  mis%PFu => cs%PFu
1081  mis%PFv => cs%PFv
1082  mis%CAu => cs%CAu
1083  mis%CAv => cs%CAv
1084  mis%pbce => cs%pbce
1085  mis%u_accel_bt => cs%u_accel_bt
1086  mis%v_accel_bt => cs%v_accel_bt
1087  mis%u_av => cs%u_av
1088  mis%v_av => cs%v_av
1089 
1090  cs%ADp => accel_diag
1091  cs%CDp => cont_diag
1092  accel_diag%diffu => cs%diffu
1093  accel_diag%diffv => cs%diffv
1094  accel_diag%PFu => cs%PFu
1095  accel_diag%PFv => cs%PFv
1096  accel_diag%CAu => cs%CAu
1097  accel_diag%CAv => cs%CAv
1098 
1099 ! Accel_diag%pbce => CS%pbce
1100 ! Accel_diag%u_accel_bt => CS%u_accel_bt ; Accel_diag%v_accel_bt => CS%v_accel_bt
1101 ! Accel_diag%u_av => CS%u_av ; Accel_diag%v_av => CS%v_av
1102 
1103  id_clock_pass_init = cpu_clock_id('(Ocean init message passing)', &
1104  grain=clock_routine)
1105 
1106  call continuity_init(time, g, gv, us, param_file, diag, cs%continuity_CSp)
1107  call coriolisadv_init(time, g, gv, us, param_file, diag, cs%ADp, cs%CoriolisAdv_CSp)
1108  if (use_tides) call tidal_forcing_init(time, g, param_file, cs%tides_CSp)
1109  call pressureforce_init(time, g, gv, us, param_file, diag, cs%PressureForce_CSp, &
1110  cs%tides_CSp)
1111  call hor_visc_init(time, g, us, param_file, diag, cs%hor_visc_CSp, meke)
1112  call vertvisc_init(mis, time, g, gv, us, param_file, diag, cs%ADp, dirs, &
1113  ntrunc, cs%vertvisc_CSp)
1114  if (.not.associated(setvisc_csp)) call mom_error(fatal, &
1115  "initialize_dyn_split_RK2 called with setVisc_CSp unassociated.")
1116  cs%set_visc_CSp => setvisc_csp
1117  call updatecfltruncationvalue(time, cs%vertvisc_CSp, &
1118  activate=is_new_run(restart_cs) )
1119 
1120  if (associated(ale_csp)) cs%ALE_CSp => ale_csp
1121  if (associated(obc)) cs%OBC => obc
1122  if (associated(update_obc_csp)) cs%update_OBC_CSp => update_obc_csp
1123 
1124  eta_rest_name = "sfc" ; if (.not.gv%Boussinesq) eta_rest_name = "p_bot"
1125  if (.not. query_initialized(cs%eta, trim(eta_rest_name), restart_cs)) then
1126  ! Estimate eta based on the layer thicknesses - h. With the Boussinesq
1127  ! approximation, eta is the free surface height anomaly, while without it
1128  ! eta is the mass of ocean per unit area. eta always has the same
1129  ! dimensions as h, either m or kg m-3.
1130  ! CS%eta(:,:) = 0.0 already from initialization.
1131  if (gv%Boussinesq) then
1132  do j=js,je ; do i=is,ie ; cs%eta(i,j) = -gv%Z_to_H * g%bathyT(i,j) ; enddo ; enddo
1133  endif
1134  do k=1,nz ; do j=js,je ; do i=is,ie
1135  cs%eta(i,j) = cs%eta(i,j) + h(i,j,k)
1136  enddo ; enddo ; enddo
1137  elseif ((gv%m_to_H_restart /= 0.0) .and. (gv%m_to_H_restart /= gv%m_to_H)) then
1138  h_rescale = gv%m_to_H / gv%m_to_H_restart
1139  do j=js,je ; do i=is,ie ; cs%eta(i,j) = h_rescale * cs%eta(i,j) ; enddo ; enddo
1140  endif
1141  ! Copy eta into an output array.
1142  do j=js,je ; do i=is,ie ; eta(i,j) = cs%eta(i,j) ; enddo ; enddo
1143 
1144  call barotropic_init(u, v, h, cs%eta, time, g, gv, us, param_file, diag, &
1145  cs%barotropic_CSp, restart_cs, calc_dtbt, cs%BT_cont, &
1146  cs%tides_CSp)
1147 
1148  if (.not. query_initialized(cs%diffu,"diffu",restart_cs) .or. &
1149  .not. query_initialized(cs%diffv,"diffv",restart_cs)) then
1150  call horizontal_viscosity(u, v, h, cs%diffu, cs%diffv, meke, varmix, &
1151  g, gv, us, cs%hor_visc_CSp, &
1152  obc=cs%OBC, bt=cs%barotropic_CSp)
1153  else
1154  if ( (us%s_to_T_restart * us%m_to_L_restart /= 0.0) .and. &
1155  (us%m_to_L * us%s_to_T_restart**2 /= us%m_to_L_restart * us%s_to_T**2) ) then
1156  accel_rescale = (us%m_to_L * us%s_to_T_restart**2) / (us%m_to_L_restart * us%s_to_T**2)
1157  do k=1,nz ; do j=js,je ; do i=g%IscB,g%IecB
1158  cs%diffu(i,j,k) = accel_rescale * cs%diffu(i,j,k)
1159  enddo ; enddo ; enddo
1160  do k=1,nz ; do j=g%JscB,g%JecB ; do i=is,ie
1161  cs%diffv(i,j,k) = accel_rescale * cs%diffv(i,j,k)
1162  enddo ; enddo ; enddo
1163  endif
1164  endif
1165 
1166  if (.not. query_initialized(cs%u_av,"u2", restart_cs) .or. &
1167  .not. query_initialized(cs%u_av,"v2", restart_cs)) then
1168  do k=1,nz ; do j=jsd,jed ; do i=isdb,iedb ; cs%u_av(i,j,k) = u(i,j,k) ; enddo ; enddo ; enddo
1169  do k=1,nz ; do j=jsdb,jedb ; do i=isd,ied ; cs%v_av(i,j,k) = v(i,j,k) ; enddo ; enddo ; enddo
1170  elseif ( (us%s_to_T_restart * us%m_to_L_restart /= 0.0) .and. &
1171  ((us%m_to_L * us%s_to_T_restart) /= (us%m_to_L_restart * us%s_to_T)) ) then
1172  vel_rescale = (us%m_to_L * us%s_to_T_restart) / (us%m_to_L_restart * us%s_to_T)
1173  do k=1,nz ; do j=jsd,jed ; do i=isdb,iedb ; cs%u_av(i,j,k) = vel_rescale * cs%u_av(i,j,k) ; enddo ; enddo ; enddo
1174  do k=1,nz ; do j=jsdb,jedb ; do i=isd,ied ; cs%v_av(i,j,k) = vel_rescale * cs%v_av(i,j,k) ; enddo ; enddo ; enddo
1175  endif
1176 
1177  ! This call is just here to initialize uh and vh.
1178  if (.not. query_initialized(uh,"uh",restart_cs) .or. &
1179  .not. query_initialized(vh,"vh",restart_cs)) then
1180  do k=1,nz ; do j=jsd,jed ; do i=isd,ied ; h_tmp(i,j,k) = h(i,j,k) ; enddo ; enddo ; enddo
1181  call continuity(u, v, h, h_tmp, uh, vh, dt, g, gv, us, cs%continuity_CSp, obc=cs%OBC)
1182  call pass_var(h_tmp, g%Domain, clock=id_clock_pass_init)
1183  do k=1,nz ; do j=jsd,jed ; do i=isd,ied
1184  cs%h_av(i,j,k) = 0.5*(h(i,j,k) + h_tmp(i,j,k))
1185  enddo ; enddo ; enddo
1186  else
1187  if (.not. query_initialized(cs%h_av,"h2",restart_cs)) then
1188  cs%h_av(:,:,:) = h(:,:,:)
1189  elseif ((gv%m_to_H_restart /= 0.0) .and. (gv%m_to_H_restart /= gv%m_to_H)) then
1190  h_rescale = gv%m_to_H / gv%m_to_H_restart
1191  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
1192  endif
1193  if ( (gv%m_to_H_restart * us%s_to_T_restart * us%m_to_L_restart /= 0.0) .and. &
1194  ((gv%m_to_H * us%m_to_L**2 * us%s_to_T_restart) /= &
1195  (gv%m_to_H_restart * us%m_to_L_restart**2 * us%s_to_T)) ) then
1196  uh_rescale = (gv%m_to_H * us%m_to_L**2 * us%s_to_T_restart) / &
1197  (gv%m_to_H_restart * us%m_to_L_restart**2 * us%s_to_T)
1198  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
1199  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
1200  endif
1201  endif
1202 
1203  call cpu_clock_begin(id_clock_pass_init)
1204  call create_group_pass(pass_av_h_uvh, cs%u_av, cs%v_av, g%Domain, halo=2)
1205  call create_group_pass(pass_av_h_uvh, cs%h_av, g%Domain, halo=2)
1206  call create_group_pass(pass_av_h_uvh, uh, vh, g%Domain, halo=2)
1207  call do_group_pass(pass_av_h_uvh, g%Domain)
1208  call cpu_clock_end(id_clock_pass_init)
1209 
1210  flux_units = get_flux_units(gv)
1211  h_convert = gv%H_to_m ; if (.not.gv%Boussinesq) h_convert = gv%H_to_kg_m2
1212  cs%id_uh = register_diag_field('ocean_model', 'uh', diag%axesCuL, time, &
1213  'Zonal Thickness Flux', flux_units, y_cell_method='sum', v_extensive=.true., &
1214  conversion=h_convert*us%L_to_m**2*us%s_to_T)
1215  cs%id_vh = register_diag_field('ocean_model', 'vh', diag%axesCvL, time, &
1216  'Meridional Thickness Flux', flux_units, x_cell_method='sum', v_extensive=.true., &
1217  conversion=h_convert*us%L_to_m**2*us%s_to_T)
1218 
1219  cs%id_CAu = register_diag_field('ocean_model', 'CAu', diag%axesCuL, time, &
1220  'Zonal Coriolis and Advective Acceleration', 'm s-2', conversion=us%L_T2_to_m_s2)
1221  cs%id_CAv = register_diag_field('ocean_model', 'CAv', diag%axesCvL, time, &
1222  'Meridional Coriolis and Advective Acceleration', 'm s-2', conversion=us%L_T2_to_m_s2)
1223  cs%id_PFu = register_diag_field('ocean_model', 'PFu', diag%axesCuL, time, &
1224  'Zonal Pressure Force Acceleration', 'm s-2', conversion=us%L_T2_to_m_s2)
1225  cs%id_PFv = register_diag_field('ocean_model', 'PFv', diag%axesCvL, time, &
1226  'Meridional Pressure Force Acceleration', 'm s-2', conversion=us%L_T2_to_m_s2)
1227 
1228  cs%id_uav = register_diag_field('ocean_model', 'uav', diag%axesCuL, time, &
1229  'Barotropic-step Averaged Zonal Velocity', 'm s-1', conversion=us%L_T_to_m_s)
1230  cs%id_vav = register_diag_field('ocean_model', 'vav', diag%axesCvL, time, &
1231  'Barotropic-step Averaged Meridional Velocity', 'm s-1', conversion=us%L_T_to_m_s)
1232 
1233  cs%id_u_BT_accel = register_diag_field('ocean_model', 'u_BT_accel', diag%axesCuL, time, &
1234  'Barotropic Anomaly Zonal Acceleration', 'm s-2', conversion=us%L_T2_to_m_s2)
1235  cs%id_v_BT_accel = register_diag_field('ocean_model', 'v_BT_accel', diag%axesCvL, time, &
1236  'Barotropic Anomaly Meridional Acceleration', 'm s-2', conversion=us%L_T2_to_m_s2)
1237 
1238  id_clock_cor = cpu_clock_id('(Ocean Coriolis & mom advection)', grain=clock_module)
1239  id_clock_continuity = cpu_clock_id('(Ocean continuity equation)', grain=clock_module)
1240  id_clock_pres = cpu_clock_id('(Ocean pressure force)', grain=clock_module)
1241  id_clock_vertvisc = cpu_clock_id('(Ocean vertical viscosity)', grain=clock_module)
1242  id_clock_horvisc = cpu_clock_id('(Ocean horizontal viscosity)', grain=clock_module)
1243  id_clock_mom_update = cpu_clock_id('(Ocean momentum increments)', grain=clock_module)
1244  id_clock_pass = cpu_clock_id('(Ocean message passing)', grain=clock_module)
1245  id_clock_btcalc = cpu_clock_id('(Ocean barotropic mode calc)', grain=clock_module)
1246  id_clock_btstep = cpu_clock_id('(Ocean barotropic mode stepping)', grain=clock_module)
1247  id_clock_btforce = cpu_clock_id('(Ocean barotropic forcing calc)', grain=clock_module)
1248 

References mom_coriolisadv::coriolisadv_init(), mom_verticalgrid::get_flux_units(), mom_hor_visc::hor_visc_init(), 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_pass_init, id_clock_pres, id_clock_vertvisc, mom_restart::is_new_run(), mom_pressureforce::pressureforce_init(), mom_tidal_forcing::tidal_forcing_init(), and mom_vert_friction::updatecfltruncationvalue().

Here is the call graph for this function:

◆ 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]hiHorizontal index structure
[in]gvocean vertical grid structure
[in]param_fileparameter file
csmodule control structure
restart_csrestart control structure
[in,out]uhzonal volume/mass transport [H L2 T-1 ~> m3 s-1 or kg s-1]
[in,out]vhmerid volume/mass transport [H L2 T-1 ~> m3 s-1 or kg s-1]

Definition at line 877 of file MOM_dynamics_split_RK2.F90.

877  type(hor_index_type), intent(in) :: HI !< Horizontal index structure
878  type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
879  type(param_file_type), intent(in) :: param_file !< parameter file
880  type(MOM_dyn_split_RK2_CS), pointer :: CS !< module control structure
881  type(MOM_restart_CS), pointer :: restart_CS !< restart control structure
882  real, dimension(SZIB_(HI),SZJ_(HI),SZK_(GV)), &
883  target, intent(inout) :: uh !< zonal volume/mass transport [H L2 T-1 ~> m3 s-1 or kg s-1]
884  real, dimension(SZI_(HI),SZJB_(HI),SZK_(GV)), &
885  target, intent(inout) :: vh !< merid volume/mass transport [H L2 T-1 ~> m3 s-1 or kg s-1]
886 
887  type(vardesc) :: vd
888  character(len=40) :: mdl = "MOM_dynamics_split_RK2" ! This module's name.
889  character(len=48) :: thickness_units, flux_units
890 
891  integer :: isd, ied, jsd, jed, nz, IsdB, IedB, JsdB, JedB
892  isd = hi%isd ; ied = hi%ied ; jsd = hi%jsd ; jed = hi%jed ; nz = gv%ke
893  isdb = hi%IsdB ; iedb = hi%IedB ; jsdb = hi%JsdB ; jedb = hi%JedB
894 
895  ! This is where a control structure specific to this module would be allocated.
896  if (associated(cs)) then
897  call mom_error(warning, "register_restarts_dyn_split_RK2 called with an associated "// &
898  "control structure.")
899  return
900  endif
901  allocate(cs)
902 
903  alloc_(cs%diffu(isdb:iedb,jsd:jed,nz)) ; cs%diffu(:,:,:) = 0.0
904  alloc_(cs%diffv(isd:ied,jsdb:jedb,nz)) ; cs%diffv(:,:,:) = 0.0
905  alloc_(cs%CAu(isdb:iedb,jsd:jed,nz)) ; cs%CAu(:,:,:) = 0.0
906  alloc_(cs%CAv(isd:ied,jsdb:jedb,nz)) ; cs%CAv(:,:,:) = 0.0
907  alloc_(cs%PFu(isdb:iedb,jsd:jed,nz)) ; cs%PFu(:,:,:) = 0.0
908  alloc_(cs%PFv(isd:ied,jsdb:jedb,nz)) ; cs%PFv(:,:,:) = 0.0
909 
910  alloc_(cs%eta(isd:ied,jsd:jed)) ; cs%eta(:,:) = 0.0
911  alloc_(cs%u_av(isdb:iedb,jsd:jed,nz)) ; cs%u_av(:,:,:) = 0.0
912  alloc_(cs%v_av(isd:ied,jsdb:jedb,nz)) ; cs%v_av(:,:,:) = 0.0
913  alloc_(cs%h_av(isd:ied,jsd:jed,nz)) ; cs%h_av(:,:,:) = gv%Angstrom_H
914 
915  thickness_units = get_thickness_units(gv)
916  flux_units = get_flux_units(gv)
917 
918  if (gv%Boussinesq) then
919  vd = var_desc("sfc",thickness_units,"Free surface Height",'h','1')
920  else
921  vd = var_desc("p_bot",thickness_units,"Bottom Pressure",'h','1')
922  endif
923  call register_restart_field(cs%eta, vd, .false., restart_cs)
924 
925  vd = var_desc("u2","m s-1","Auxiliary Zonal velocity",'u','L')
926  call register_restart_field(cs%u_av, vd, .false., restart_cs)
927 
928  vd = var_desc("v2","m s-1","Auxiliary Meridional velocity",'v','L')
929  call register_restart_field(cs%v_av, vd, .false., restart_cs)
930 
931  vd = var_desc("h2",thickness_units,"Auxiliary Layer Thickness",'h','L')
932  call register_restart_field(cs%h_av, vd, .false., restart_cs)
933 
934  vd = var_desc("uh",flux_units,"Zonal thickness flux",'u','L')
935  call register_restart_field(uh, vd, .false., restart_cs)
936 
937  vd = var_desc("vh",flux_units,"Meridional thickness flux",'v','L')
938  call register_restart_field(vh, vd, .false., restart_cs)
939 
940  vd = var_desc("diffu","m s-2","Zonal horizontal viscous acceleration",'u','L')
941  call register_restart_field(cs%diffu, vd, .false., restart_cs)
942 
943  vd = var_desc("diffv","m s-2","Meridional horizontal viscous acceleration",'v','L')
944  call register_restart_field(cs%diffv, vd, .false., restart_cs)
945 
946  call register_barotropic_restarts(hi, gv, param_file, cs%barotropic_CSp, &
947  restart_cs)
948 

References mom_verticalgrid::get_flux_units(), mom_barotropic::register_barotropic_restarts(), and mom_io::var_desc().

Here is the call graph for this function:

◆ 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]gocean grid structure
[in]gvocean vertical grid structure
[in]usA dimensional unit scaling type
[in,out]uzonal velocity [L T-1 ~> m s-1]
[in,out]vmerid velocity [L T-1 ~> m s-1]
[in,out]hlayer thickness [H ~> m or kg m-2]
[in]tvthermodynamic type
[in,out]viscvertical visc, bottom drag, and related
[in]time_localmodel time at end of time step
[in]dttime step [T ~> s]
[in]forcesA structure with the driving mechanical forces
p_surf_beginsurf pressure at start of this dynamic time step [Pa]
p_surf_endsurf pressure at end of this dynamic time step [Pa]
[in,out]uhzonal volume/mass transport
[in,out]vhmerid volume/mass transport
[in,out]uhtraccumulatated zonal volume/mass transport
[in,out]vhtraccumulatated merid volume/mass transport
[out]eta_avfree surface height or column mass time averaged over time step [H ~> m or kg m-2]
csmodule control structure
[in]calc_dtbtif true, recalculate barotropic time step
varmixspecify the spatially varying viscosities
mekerelated to mesoscale eddy kinetic energy param
thickness_diffuse_cspPointer to a structure containing interface height diffusivities
wavesA pointer to a structure containing fields related to the surface wave conditions

Definition at line 239 of file MOM_dynamics_split_RK2.F90.

239  type(ocean_grid_type), intent(inout) :: G !< ocean grid structure
240  type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
241  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
242  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
243  target, intent(inout) :: u !< zonal velocity [L T-1 ~> m s-1]
244  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
245  target, intent(inout) :: v !< merid velocity [L T-1 ~> m s-1]
246  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
247  intent(inout) :: h !< layer thickness [H ~> m or kg m-2]
248  type(thermo_var_ptrs), intent(in) :: tv !< thermodynamic type
249  type(vertvisc_type), intent(inout) :: visc !< vertical visc, bottom drag, and related
250  type(time_type), intent(in) :: Time_local !< model time at end of time step
251  real, intent(in) :: dt !< time step [T ~> s]
252  type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces
253  real, dimension(:,:), pointer :: p_surf_begin !< surf pressure at start of this dynamic
254  !! time step [Pa]
255  real, dimension(:,:), pointer :: p_surf_end !< surf pressure at end of this dynamic
256  !! time step [Pa]
257  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
258  target, intent(inout) :: uh !< zonal volume/mass transport
259  !! [H L2 T-1 ~> m3 s-1 or kg s-1]
260  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
261  target, intent(inout) :: vh !< merid volume/mass transport
262  !! [H L2 T-1 ~> m3 s-1 or kg s-1]
263  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
264  intent(inout) :: uhtr !< accumulatated zonal volume/mass transport
265  !! since last tracer advection [H L2 ~> m3 or kg]
266  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
267  intent(inout) :: vhtr !< accumulatated merid volume/mass transport
268  !! since last tracer advection [H L2 ~> m3 or kg]
269  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: eta_av !< free surface height or column mass time
270  !! averaged over time step [H ~> m or kg m-2]
271  type(MOM_dyn_split_RK2_CS), pointer :: CS !< module control structure
272  logical, intent(in) :: calc_dtbt !< if true, recalculate barotropic time step
273  type(VarMix_CS), pointer :: VarMix !< specify the spatially varying viscosities
274  type(MEKE_type), pointer :: MEKE !< related to mesoscale eddy kinetic energy param
275  type(thickness_diffuse_CS), pointer :: thickness_diffuse_CSp!< Pointer to a structure containing
276  !! interface height diffusivities
277  type(wave_parameters_CS), optional, pointer :: Waves !< A pointer to a structure containing
278  !! fields related to the surface wave conditions
279 
280  ! local variables
281  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: up ! Predicted zonal velocity [L T-1 ~> m s-1].
282  real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: vp ! Predicted meridional velocity [L T-1 ~> m s-1].
283  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: hp ! Predicted thickness [H ~> m or kg m-2].
284 
285  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: u_bc_accel
286  real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: v_bc_accel
287  ! u_bc_accel and v_bc_accel are the summed baroclinic accelerations of each
288  ! layer calculated by the non-barotropic part of the model [L T-2 ~> m s-2].
289 
290  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), target :: uh_in
291  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), target :: vh_in
292  ! uh_in and vh_in are the zonal or meridional mass transports that would be
293  ! obtained using the initial velocities [H L2 T-1 ~> m3 s-1 or kg s-1].
294 
295  real, dimension(SZIB_(G),SZJ_(G)) :: uhbt_out
296  real, dimension(SZI_(G),SZJB_(G)) :: vhbt_out
297  ! uhbt_out and vhbt_out are the vertically summed transports from the
298  ! barotropic solver based on its final velocities [H m2 s-1 ~> m3 s-1 or kg s-1].
299 
300  real, dimension(SZI_(G),SZJ_(G)) :: eta_pred
301  ! eta_pred is the predictor value of the free surface height or column mass,
302  ! [H ~> m or kg m-2].
303 
304  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: u_old_rad_OBC
305  real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: v_old_rad_OBC
306  ! u_old_rad_OBC and v_old_rad_OBC are the starting velocities, which are
307  ! saved for use in the Flather open boundary condition code [L T-1 ~> m s-1].
308 
309  real :: Pa_to_eta ! A factor that converts pressures to the units of eta.
310  real, pointer, dimension(:,:) :: &
311  p_surf => null(), eta_pf_start => null(), &
312  taux_bot => null(), tauy_bot => null(), &
313  eta => null()
314 
315  real, pointer, dimension(:,:,:) :: &
316  uh_ptr => null(), u_ptr => null(), vh_ptr => null(), v_ptr => null(), &
317  u_av, & ! The zonal velocity time-averaged over a time step [L T-1 ~> m s-1].
318  v_av, & ! The meridional velocity time-averaged over a time step [L T-1 ~> m s-1].
319  h_av ! The layer thickness time-averaged over a time step [H ~> m or kg m-2].
320  real :: dt_pred ! The time step for the predictor part of the baroclinic time stepping [T ~> s].
321 
322  logical :: dyn_p_surf
323  logical :: BT_cont_BT_thick ! If true, use the BT_cont_type to estimate the
324  ! relative weightings of the layers in calculating
325  ! the barotropic accelerations.
326  !---For group halo pass
327  logical :: showCallTree, sym
328 
329  integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
330  integer :: cont_stencil
331  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
332  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
333  u_av => cs%u_av ; v_av => cs%v_av ; h_av => cs%h_av ; eta => cs%eta
334 
335  sym=.false.;if (g%Domain%symmetric) sym=.true. ! switch to include symmetric domain in checksums
336 
337  showcalltree = calltree_showquery()
338  if (showcalltree) call calltree_enter("step_MOM_dyn_split_RK2(), MOM_dynamics_split_RK2.F90")
339 
340  !$OMP parallel do default(shared)
341  do k=1,nz
342  do j=g%jsd,g%jed ; do i=g%isdB,g%iedB ; up(i,j,k) = 0.0 ; enddo ; enddo
343  do j=g%jsdB,g%jedB ; do i=g%isd,g%ied ; vp(i,j,k) = 0.0 ; enddo ; enddo
344  do j=g%jsd,g%jed ; do i=g%isd,g%ied ; hp(i,j,k) = h(i,j,k) ; enddo ; enddo
345  enddo
346 
347  ! Update CFL truncation value as function of time
348  call updatecfltruncationvalue(time_local, cs%vertvisc_CSp)
349 
350  if (cs%debug) then
351  call mom_state_chksum("Start predictor ", u, v, h, uh, vh, g, gv, us, symmetric=sym)
352  call check_redundant("Start predictor u ", u, v, g)
353  call check_redundant("Start predictor uh ", uh, vh, g)
354  endif
355 
356  dyn_p_surf = associated(p_surf_begin) .and. associated(p_surf_end)
357  if (dyn_p_surf) then
358  p_surf => p_surf_end
359  call safe_alloc_ptr(eta_pf_start,g%isd,g%ied,g%jsd,g%jed)
360  eta_pf_start(:,:) = 0.0
361  else
362  p_surf => forces%p_surf
363  endif
364 
365  if (associated(cs%OBC)) then
366  if (cs%debug_OBC) call open_boundary_test_extern_h(g, gv, cs%OBC, h)
367 
368  do k=1,nz ; do j=g%jsd,g%jed ; do i=g%IsdB,g%IedB
369  u_old_rad_obc(i,j,k) = u_av(i,j,k)
370  enddo ; enddo ; enddo
371  do k=1,nz ; do j=g%JsdB,g%JedB ; do i=g%isd,g%ied
372  v_old_rad_obc(i,j,k) = v_av(i,j,k)
373  enddo ; enddo ; enddo
374  endif
375 
376  bt_cont_bt_thick = .false.
377  if (associated(cs%BT_cont)) bt_cont_bt_thick = &
378  (allocated(cs%BT_cont%h_u) .and. allocated(cs%BT_cont%h_v))
379 
380  if (cs%split_bottom_stress) then
381  taux_bot => cs%taux_bot ; tauy_bot => cs%tauy_bot
382  endif
383 
384  !--- begin set up for group halo pass
385 
386  cont_stencil = continuity_stencil(cs%continuity_CSp)
387  !### Apart from circle_OBCs halo for eta could be 1, but halo>=3 is required
388  !### to match circle_OBCs solutions. Why?
389  call cpu_clock_begin(id_clock_pass)
390  call create_group_pass(cs%pass_eta, eta, g%Domain) !### , halo=1)
391  call create_group_pass(cs%pass_visc_rem, cs%visc_rem_u, cs%visc_rem_v, g%Domain, &
392  to_all+scalar_pair, cgrid_ne, halo=max(1,cont_stencil))
393  call create_group_pass(cs%pass_uvp, up, vp, g%Domain, halo=max(1,cont_stencil))
394  call create_group_pass(cs%pass_hp_uv, hp, g%Domain, halo=2)
395  call create_group_pass(cs%pass_hp_uv, u_av, v_av, g%Domain, halo=2)
396  call create_group_pass(cs%pass_hp_uv, uh(:,:,:), vh(:,:,:), g%Domain, halo=2)
397 
398  call create_group_pass(cs%pass_uv, u, v, g%Domain, halo=max(2,cont_stencil))
399  call create_group_pass(cs%pass_h, h, g%Domain, halo=max(2,cont_stencil))
400  call create_group_pass(cs%pass_av_uvh, u_av, v_av, g%Domain, halo=2)
401  call create_group_pass(cs%pass_av_uvh, uh(:,:,:), vh(:,:,:), g%Domain, halo=2)
402  call cpu_clock_end(id_clock_pass)
403  !--- end set up for group halo pass
404 
405 
406 ! PFu = d/dx M(h,T,S)
407 ! pbce = dM/deta
408  if (cs%begw == 0.0) call enable_averages(dt, time_local, cs%diag)
409  call cpu_clock_begin(id_clock_pres)
410  call pressureforce(h, tv, cs%PFu, cs%PFv, g, gv, us, cs%PressureForce_CSp, &
411  cs%ALE_CSp, p_surf, cs%pbce, cs%eta_PF)
412  if (dyn_p_surf) then
413  pa_to_eta = 1.0 / gv%H_to_Pa
414  !$OMP parallel do default(shared)
415  do j=jsq,jeq+1 ; do i=isq,ieq+1
416  eta_pf_start(i,j) = cs%eta_PF(i,j) - pa_to_eta * &
417  (p_surf_begin(i,j) - p_surf_end(i,j))
418  enddo ; enddo
419  endif
420  call cpu_clock_end(id_clock_pres)
421  call disable_averaging(cs%diag)
422  if (showcalltree) call calltree_waypoint("done with PressureForce (step_MOM_dyn_split_RK2)")
423 
424  if (associated(cs%OBC)) then; if (cs%OBC%update_OBC) then
425  call update_obc_data(cs%OBC, g, gv, us, tv, h, cs%update_OBC_CSp, time_local)
426  endif; endif
427  if (associated(cs%OBC) .and. cs%debug_OBC) &
428  call open_boundary_zero_normal_flow(cs%OBC, g, cs%PFu, cs%PFv)
429 
430  if (g%nonblocking_updates) &
431  call start_group_pass(cs%pass_eta, g%Domain, clock=id_clock_pass)
432 
433 ! CAu = -(f+zeta_av)/h_av vh + d/dx KE_av
434  call cpu_clock_begin(id_clock_cor)
435  call coradcalc(u_av, v_av, h_av, uh, vh, cs%CAu, cs%CAv, cs%OBC, cs%ADp, &
436  g, gv, us, cs%CoriolisAdv_CSp)
437  call cpu_clock_end(id_clock_cor)
438  if (showcalltree) call calltree_waypoint("done with CorAdCalc (step_MOM_dyn_split_RK2)")
439 
440 ! u_bc_accel = CAu + PFu + diffu(u[n-1])
441  call cpu_clock_begin(id_clock_btforce)
442  !$OMP parallel do default(shared)
443  do k=1,nz
444  do j=js,je ; do i=isq,ieq
445  u_bc_accel(i,j,k) = (cs%CAu(i,j,k) + cs%PFu(i,j,k)) + cs%diffu(i,j,k)
446  enddo ; enddo
447  do j=jsq,jeq ; do i=is,ie
448  v_bc_accel(i,j,k) = (cs%CAv(i,j,k) + cs%PFv(i,j,k)) + cs%diffv(i,j,k)
449  enddo ; enddo
450  enddo
451  if (associated(cs%OBC)) then
452  call open_boundary_zero_normal_flow(cs%OBC, g, u_bc_accel, v_bc_accel)
453  endif
454  call cpu_clock_end(id_clock_btforce)
455 
456  if (cs%debug) then
457  call mom_accel_chksum("pre-btstep accel", cs%CAu, cs%CAv, cs%PFu, cs%PFv, &
458  cs%diffu, cs%diffv, g, gv, us, cs%pbce, u_bc_accel, v_bc_accel, &
459  symmetric=sym)
460  call check_redundant("pre-btstep CS%Ca ", cs%Cau, cs%Cav, g)
461  call check_redundant("pre-btstep CS%PF ", cs%PFu, cs%PFv, g)
462  call check_redundant("pre-btstep CS%diff ", cs%diffu, cs%diffv, g)
463  call check_redundant("pre-btstep u_bc_accel ", u_bc_accel, v_bc_accel, g)
464  endif
465 
466  call cpu_clock_begin(id_clock_vertvisc)
467  !$OMP parallel do default(shared)
468  do k=1,nz
469  do j=js,je ; do i=isq,ieq
470  up(i,j,k) = g%mask2dCu(i,j) * (u(i,j,k) + dt * u_bc_accel(i,j,k))
471  enddo ; enddo
472  do j=jsq,jeq ; do i=is,ie
473  vp(i,j,k) = g%mask2dCv(i,j) * (v(i,j,k) + dt * v_bc_accel(i,j,k))
474  enddo ; enddo
475  enddo
476 
477  call enable_averages(dt, time_local, cs%diag)
478  call set_viscous_ml(u, v, h, tv, forces, visc, dt, g, gv, us, &
479  cs%set_visc_CSp)
480  call disable_averaging(cs%diag)
481 
482  if (cs%debug) then
483  call uvchksum("before vertvisc: up", up, vp, g%HI, haloshift=0, symmetric=sym, scale=us%L_T_to_m_s)
484  endif
485  call vertvisc_coef(up, vp, h, forces, visc, dt, g, gv, us, cs%vertvisc_CSp, cs%OBC)
486  call vertvisc_remnant(visc, cs%visc_rem_u, cs%visc_rem_v, dt, g, gv, us, cs%vertvisc_CSp)
487  call cpu_clock_end(id_clock_vertvisc)
488  if (showcalltree) call calltree_waypoint("done with vertvisc_coef (step_MOM_dyn_split_RK2)")
489 
490 
491  call cpu_clock_begin(id_clock_pass)
492  if (g%nonblocking_updates) then
493  call complete_group_pass(cs%pass_eta, g%Domain)
494  call start_group_pass(cs%pass_visc_rem, g%Domain)
495  else
496  call do_group_pass(cs%pass_eta, g%Domain)
497  call do_group_pass(cs%pass_visc_rem, g%Domain)
498  endif
499  call cpu_clock_end(id_clock_pass)
500 
501  call cpu_clock_begin(id_clock_btcalc)
502  ! Calculate the relative layer weights for determining barotropic quantities.
503  if (.not.bt_cont_bt_thick) &
504  call btcalc(h, g, gv, cs%barotropic_CSp, obc=cs%OBC)
505  call bt_mass_source(h, eta, .true., g, gv, cs%barotropic_CSp)
506  call cpu_clock_end(id_clock_btcalc)
507 
508  if (g%nonblocking_updates) &
509  call complete_group_pass(cs%pass_visc_rem, g%Domain, clock=id_clock_pass)
510 
511 ! u_accel_bt = layer accelerations due to barotropic solver
512  if (associated(cs%BT_cont) .or. cs%BT_use_layer_fluxes) then
513  call cpu_clock_begin(id_clock_continuity)
514  call continuity(u, v, h, hp, uh_in, vh_in, dt, g, gv, us, cs%continuity_CSp, &
515  obc=cs%OBC, visc_rem_u=cs%visc_rem_u, visc_rem_v=cs%visc_rem_v, bt_cont=cs%BT_cont)
516  call cpu_clock_end(id_clock_continuity)
517  if (bt_cont_bt_thick) then
518  call btcalc(h, g, gv, cs%barotropic_CSp, cs%BT_cont%h_u, cs%BT_cont%h_v, &
519  obc=cs%OBC)
520  endif
521  if (showcalltree) call calltree_waypoint("done with continuity[BT_cont] (step_MOM_dyn_split_RK2)")
522  endif
523 
524  if (cs%BT_use_layer_fluxes) then
525  uh_ptr => uh_in ; vh_ptr => vh_in; u_ptr => u ; v_ptr => v
526  endif
527 
528  call cpu_clock_begin(id_clock_btstep)
529  if (calc_dtbt) call set_dtbt(g, gv, us, cs%barotropic_CSp, eta, cs%pbce)
530  if (showcalltree) call calltree_enter("btstep(), MOM_barotropic.F90")
531  ! This is the predictor step call to btstep.
532  call btstep(u, v, eta, dt, u_bc_accel, v_bc_accel, forces, cs%pbce, cs%eta_PF, &
533  u_av, v_av, cs%u_accel_bt, cs%v_accel_bt, eta_pred, cs%uhbt, cs%vhbt, &
534  g, gv, us, cs%barotropic_CSp, cs%visc_rem_u, cs%visc_rem_v, &
535  obc=cs%OBC, bt_cont=cs%BT_cont, eta_pf_start=eta_pf_start, &
536  taux_bot=taux_bot, tauy_bot=tauy_bot, &
537  uh0=uh_ptr, vh0=vh_ptr, u_uh0=u_ptr, v_vh0=v_ptr)
538  if (showcalltree) call calltree_leave("btstep()")
539  call cpu_clock_end(id_clock_btstep)
540 
541 ! up = u + dt_pred*( u_bc_accel + u_accel_bt )
542  dt_pred = dt * cs%be
543  call cpu_clock_begin(id_clock_mom_update)
544 
545  !$OMP parallel do default(shared)
546  do k=1,nz
547  do j=jsq,jeq ; do i=is,ie
548  vp(i,j,k) = g%mask2dCv(i,j) * (v(i,j,k) + dt_pred * &
549  (v_bc_accel(i,j,k) + cs%v_accel_bt(i,j,k)))
550  enddo ; enddo
551  do j=js,je ; do i=isq,ieq
552  up(i,j,k) = g%mask2dCu(i,j) * (u(i,j,k) + dt_pred * &
553  (u_bc_accel(i,j,k) + cs%u_accel_bt(i,j,k)))
554  enddo ; enddo
555  enddo
556  call cpu_clock_end(id_clock_mom_update)
557 
558  if (cs%debug) then
559  call uvchksum("Predictor 1 [uv]", up, vp, g%HI, haloshift=0, symmetric=sym, scale=us%L_T_to_m_s)
560  call hchksum(h, "Predictor 1 h", g%HI, haloshift=1, scale=gv%H_to_m)
561  call uvchksum("Predictor 1 [uv]h", uh, vh, g%HI,haloshift=2, &
562  symmetric=sym, scale=gv%H_to_m*us%L_to_m**2*us%s_to_T)
563 ! call MOM_state_chksum("Predictor 1", up, vp, h, uh, vh, G, GV, US, haloshift=1)
564  call mom_accel_chksum("Predictor accel", cs%CAu, cs%CAv, cs%PFu, cs%PFv, &
565  cs%diffu, cs%diffv, g, gv, us, cs%pbce, cs%u_accel_bt, cs%v_accel_bt, symmetric=sym)
566  call mom_state_chksum("Predictor 1 init", u, v, h, uh, vh, g, gv, us, haloshift=2, &
567  symmetric=sym)
568  call check_redundant("Predictor 1 up", up, vp, g)
569  call check_redundant("Predictor 1 uh", uh, vh, g)
570  endif
571 
572 ! up <- up + dt_pred d/dz visc d/dz up
573 ! u_av <- u_av + dt_pred d/dz visc d/dz u_av
574  call cpu_clock_begin(id_clock_vertvisc)
575  if (cs%debug) then
576  call uvchksum("0 before vertvisc: [uv]p", up, vp, g%HI,haloshift=0, symmetric=sym, scale=us%L_T_to_m_s)
577  endif
578  call vertvisc_coef(up, vp, h, forces, visc, dt_pred, g, gv, us, cs%vertvisc_CSp, &
579  cs%OBC)
580  call vertvisc(up, vp, h, forces, visc, dt_pred, cs%OBC, cs%ADp, cs%CDp, g, &
581  gv, us, cs%vertvisc_CSp, cs%taux_bot, cs%tauy_bot, waves=waves)
582  if (showcalltree) call calltree_waypoint("done with vertvisc (step_MOM_dyn_split_RK2)")
583  if (g%nonblocking_updates) then
584  call cpu_clock_end(id_clock_vertvisc)
585  call start_group_pass(cs%pass_uvp, g%Domain, clock=id_clock_pass)
586  call cpu_clock_begin(id_clock_vertvisc)
587  endif
588  call vertvisc_remnant(visc, cs%visc_rem_u, cs%visc_rem_v, dt_pred, g, gv, us, cs%vertvisc_CSp)
589  call cpu_clock_end(id_clock_vertvisc)
590 
591  call do_group_pass(cs%pass_visc_rem, g%Domain, clock=id_clock_pass)
592  if (g%nonblocking_updates) then
593  call complete_group_pass(cs%pass_uvp, g%Domain, clock=id_clock_pass)
594  else
595  call do_group_pass(cs%pass_uvp, g%Domain, clock=id_clock_pass)
596  endif
597 
598  ! uh = u_av * h
599  ! hp = h + dt * div . uh
600  call cpu_clock_begin(id_clock_continuity)
601  call continuity(up, vp, h, hp, uh, vh, dt, g, gv, us, cs%continuity_CSp, &
602  cs%uhbt, cs%vhbt, cs%OBC, cs%visc_rem_u, cs%visc_rem_v, &
603  u_av, v_av, bt_cont=cs%BT_cont)
604  call cpu_clock_end(id_clock_continuity)
605  if (showcalltree) call calltree_waypoint("done with continuity (step_MOM_dyn_split_RK2)")
606 
607  call do_group_pass(cs%pass_hp_uv, g%Domain, clock=id_clock_pass)
608 
609  if (associated(cs%OBC)) then
610 
611  if (cs%debug) &
612  call uvchksum("Pre OBC avg [uv]", u_av, v_av, g%HI, haloshift=1, symmetric=sym, scale=us%L_T_to_m_s)
613 
614  call radiation_open_bdry_conds(cs%OBC, u_av, u_old_rad_obc, v_av, v_old_rad_obc, g, us, dt_pred)
615 
616  if (cs%debug) &
617  call uvchksum("Post OBC avg [uv]", u_av, v_av, g%HI, haloshift=1, symmetric=sym, scale=us%L_T_to_m_s)
618 
619  ! These should be done with a pass that excludes uh & vh.
620 ! call do_group_pass(CS%pass_hp_uv, G%Domain, clock=id_clock_pass)
621  endif
622 
623  if (g%nonblocking_updates) then
624  call start_group_pass(cs%pass_av_uvh, g%Domain, clock=id_clock_pass)
625  endif
626 
627  ! h_av = (h + hp)/2
628  !$OMP parallel do default(shared)
629  do k=1,nz ; do j=js-2,je+2 ; do i=is-2,ie+2
630  h_av(i,j,k) = 0.5*(h(i,j,k) + hp(i,j,k))
631  enddo ; enddo ; enddo
632 
633  ! The correction phase of the time step starts here.
634  call enable_averages(dt, time_local, cs%diag)
635 
636  ! Calculate a revised estimate of the free-surface height correction to be
637  ! used in the next call to btstep. This call is at this point so that
638  ! hp can be changed if CS%begw /= 0.
639  ! eta_cor = ... (hidden inside CS%barotropic_CSp)
640  call cpu_clock_begin(id_clock_btcalc)
641  call bt_mass_source(hp, eta_pred, .false., g, gv, cs%barotropic_CSp)
642  call cpu_clock_end(id_clock_btcalc)
643 
644  if (cs%begw /= 0.0) then
645  ! hp <- (1-begw)*h_in + begw*hp
646  ! Back up hp to the value it would have had after a time-step of
647  ! begw*dt. hp is not used again until recalculated by continuity.
648  !$OMP parallel do default(shared)
649  do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1
650  hp(i,j,k) = (1.0-cs%begw)*h(i,j,k) + cs%begw*hp(i,j,k)
651  enddo ; enddo ; enddo
652 
653  ! PFu = d/dx M(hp,T,S)
654  ! pbce = dM/deta
655  call cpu_clock_begin(id_clock_pres)
656  call pressureforce(hp, tv, cs%PFu, cs%PFv, g, gv, us, cs%PressureForce_CSp, &
657  cs%ALE_CSp, p_surf, cs%pbce, cs%eta_PF)
658  call cpu_clock_end(id_clock_pres)
659  if (showcalltree) call calltree_waypoint("done with PressureForce[hp=(1-b).h+b.h] (step_MOM_dyn_split_RK2)")
660  endif
661 
662  if (g%nonblocking_updates) &
663  call complete_group_pass(cs%pass_av_uvh, g%Domain, clock=id_clock_pass)
664 
665  if (bt_cont_bt_thick) then
666  call btcalc(h, g, gv, cs%barotropic_CSp, cs%BT_cont%h_u, cs%BT_cont%h_v, &
667  obc=cs%OBC)
668  if (showcalltree) call calltree_waypoint("done with btcalc[BT_cont_BT_thick] (step_MOM_dyn_split_RK2)")
669  endif
670 
671  if (cs%debug) then
672  call mom_state_chksum("Predictor ", up, vp, hp, uh, vh, g, gv, us, symmetric=sym)
673  call uvchksum("Predictor avg [uv]", u_av, v_av, g%HI, haloshift=1, symmetric=sym, scale=us%L_T_to_m_s)
674  call hchksum(h_av, "Predictor avg h", g%HI, haloshift=0, scale=gv%H_to_m)
675  ! call MOM_state_chksum("Predictor avg ", u_av, v_av, h_av, uh, vh, G, GV, US)
676  call check_redundant("Predictor up ", up, vp, g)
677  call check_redundant("Predictor uh ", uh, vh, g)
678  endif
679 
680 ! diffu = horizontal viscosity terms (u_av)
681  call cpu_clock_begin(id_clock_horvisc)
682  call horizontal_viscosity(u_av, v_av, h_av, cs%diffu, cs%diffv, &
683  meke, varmix, g, gv, us, cs%hor_visc_CSp, &
684  obc=cs%OBC, bt=cs%barotropic_CSp)
685  call cpu_clock_end(id_clock_horvisc)
686  if (showcalltree) call calltree_waypoint("done with horizontal_viscosity (step_MOM_dyn_split_RK2)")
687 
688 ! CAu = -(f+zeta_av)/h_av vh + d/dx KE_av
689  call cpu_clock_begin(id_clock_cor)
690  call coradcalc(u_av, v_av, h_av, uh, vh, cs%CAu, cs%CAv, cs%OBC, cs%ADp, &
691  g, gv, us, cs%CoriolisAdv_CSp)
692  call cpu_clock_end(id_clock_cor)
693  if (showcalltree) call calltree_waypoint("done with CorAdCalc (step_MOM_dyn_split_RK2)")
694 
695 ! Calculate the momentum forcing terms for the barotropic equations.
696 
697 ! u_bc_accel = CAu + PFu + diffu(u[n-1])
698  call cpu_clock_begin(id_clock_btforce)
699  !$OMP parallel do default(shared)
700  do k=1,nz
701  do j=js,je ; do i=isq,ieq
702  u_bc_accel(i,j,k) = (cs%Cau(i,j,k) + cs%PFu(i,j,k)) + cs%diffu(i,j,k)
703  enddo ; enddo
704  do j=jsq,jeq ; do i=is,ie
705  v_bc_accel(i,j,k) = (cs%Cav(i,j,k) + cs%PFv(i,j,k)) + cs%diffv(i,j,k)
706  enddo ; enddo
707  enddo
708  if (associated(cs%OBC)) then
709  call open_boundary_zero_normal_flow(cs%OBC, g, u_bc_accel, v_bc_accel)
710  endif
711  call cpu_clock_end(id_clock_btforce)
712 
713  if (cs%debug) then
714  call mom_accel_chksum("corr pre-btstep accel", cs%CAu, cs%CAv, cs%PFu, cs%PFv, &
715  cs%diffu, cs%diffv, g, gv, us, cs%pbce, u_bc_accel, v_bc_accel, &
716  symmetric=sym)
717  call check_redundant("corr pre-btstep CS%Ca ", cs%Cau, cs%Cav, g)
718  call check_redundant("corr pre-btstep CS%PF ", cs%PFu, cs%PFv, g)
719  call check_redundant("corr pre-btstep CS%diff ", cs%diffu, cs%diffv, g)
720  call check_redundant("corr pre-btstep u_bc_accel ", u_bc_accel, v_bc_accel, g)
721  endif
722 
723  ! u_accel_bt = layer accelerations due to barotropic solver
724  ! pbce = dM/deta
725  call cpu_clock_begin(id_clock_btstep)
726  if (cs%BT_use_layer_fluxes) then
727  uh_ptr => uh ; vh_ptr => vh ; u_ptr => u_av ; v_ptr => v_av
728  endif
729 
730  if (showcalltree) call calltree_enter("btstep(), MOM_barotropic.F90")
731  ! This is the corrector step call to btstep.
732  call btstep(u, v, eta, dt, u_bc_accel, v_bc_accel, forces, cs%pbce, &
733  cs%eta_PF, u_av, v_av, cs%u_accel_bt, cs%v_accel_bt, &
734  eta_pred, cs%uhbt, cs%vhbt, g, gv, us, cs%barotropic_CSp, &
735  cs%visc_rem_u, cs%visc_rem_v, etaav=eta_av, obc=cs%OBC, &
736  bt_cont = cs%BT_cont, eta_pf_start=eta_pf_start, &
737  taux_bot=taux_bot, tauy_bot=tauy_bot, &
738  uh0=uh_ptr, vh0=vh_ptr, u_uh0=u_ptr, v_vh0=v_ptr)
739  do j=js,je ; do i=is,ie ; eta(i,j) = eta_pred(i,j) ; enddo ; enddo
740 
741  call cpu_clock_end(id_clock_btstep)
742  if (showcalltree) call calltree_leave("btstep()")
743 
744  if (cs%debug) then
745  call check_redundant("u_accel_bt ", cs%u_accel_bt, cs%v_accel_bt, g)
746  endif
747 
748  ! u = u + dt*( u_bc_accel + u_accel_bt )
749  call cpu_clock_begin(id_clock_mom_update)
750  !$OMP parallel do default(shared)
751  do k=1,nz
752  do j=js,je ; do i=isq,ieq
753  u(i,j,k) = g%mask2dCu(i,j) * (u(i,j,k) + dt * &
754  (u_bc_accel(i,j,k) + cs%u_accel_bt(i,j,k)))
755  enddo ; enddo
756  do j=jsq,jeq ; do i=is,ie
757  v(i,j,k) = g%mask2dCv(i,j) * (v(i,j,k) + dt * &
758  (v_bc_accel(i,j,k) + cs%v_accel_bt(i,j,k)))
759  enddo ; enddo
760  enddo
761  call cpu_clock_end(id_clock_mom_update)
762 
763  if (cs%debug) then
764  call uvchksum("Corrector 1 [uv]", u, v, g%HI,haloshift=0, symmetric=sym, scale=us%L_T_to_m_s)
765  call hchksum(h, "Corrector 1 h", g%HI, haloshift=2, scale=gv%H_to_m)
766  call uvchksum("Corrector 1 [uv]h", uh, vh, g%HI, haloshift=2, &
767  symmetric=sym, scale=gv%H_to_m*us%L_to_m**2*us%s_to_T)
768  ! call MOM_state_chksum("Corrector 1", u, v, h, uh, vh, G, GV, US, haloshift=1)
769  call mom_accel_chksum("Corrector accel", cs%CAu, cs%CAv, cs%PFu, cs%PFv, &
770  cs%diffu, cs%diffv, g, gv, us, cs%pbce, cs%u_accel_bt, cs%v_accel_bt, &
771  symmetric=sym)
772  endif
773 
774  ! u <- u + dt d/dz visc d/dz u
775  ! u_av <- u_av + dt d/dz visc d/dz u_av
776  call cpu_clock_begin(id_clock_vertvisc)
777  call vertvisc_coef(u, v, h, forces, visc, dt, g, gv, us, cs%vertvisc_CSp, cs%OBC)
778  call vertvisc(u, v, h, forces, visc, dt, cs%OBC, cs%ADp, cs%CDp, g, gv, us, &
779  cs%vertvisc_CSp, cs%taux_bot, cs%tauy_bot,waves=waves)
780  if (g%nonblocking_updates) then
781  call cpu_clock_end(id_clock_vertvisc)
782  call start_group_pass(cs%pass_uv, g%Domain, clock=id_clock_pass)
783  call cpu_clock_begin(id_clock_vertvisc)
784  endif
785  call vertvisc_remnant(visc, cs%visc_rem_u, cs%visc_rem_v, dt, g, gv, us, cs%vertvisc_CSp)
786  call cpu_clock_end(id_clock_vertvisc)
787  if (showcalltree) call calltree_waypoint("done with vertvisc (step_MOM_dyn_split_RK2)")
788 
789 ! Later, h_av = (h_in + h_out)/2, but for now use h_av to store h_in.
790  !$OMP parallel do default(shared)
791  do k=1,nz ; do j=js-2,je+2 ; do i=is-2,ie+2
792  h_av(i,j,k) = h(i,j,k)
793  enddo ; enddo ; enddo
794 
795  call do_group_pass(cs%pass_visc_rem, g%Domain, clock=id_clock_pass)
796  if (g%nonblocking_updates) then
797  call complete_group_pass(cs%pass_uv, g%Domain, clock=id_clock_pass)
798  else
799  call do_group_pass(cs%pass_uv, g%Domain, clock=id_clock_pass)
800  endif
801 
802  ! uh = u_av * h
803  ! h = h + dt * div . uh
804  ! u_av and v_av adjusted so their mass transports match uhbt and vhbt.
805  call cpu_clock_begin(id_clock_continuity)
806  call continuity(u, v, h, h, uh, vh, dt, g, gv, us, cs%continuity_CSp, &
807  cs%uhbt, cs%vhbt, cs%OBC, cs%visc_rem_u, cs%visc_rem_v, u_av, v_av)
808  call cpu_clock_end(id_clock_continuity)
809  call do_group_pass(cs%pass_h, g%Domain, clock=id_clock_pass)
810  ! Whenever thickness changes let the diag manager know, target grids
811  ! for vertical remapping may need to be regenerated.
812  call diag_update_remap_grids(cs%diag)
813  if (showcalltree) call calltree_waypoint("done with continuity (step_MOM_dyn_split_RK2)")
814 
815  if (g%nonblocking_updates) then
816  call start_group_pass(cs%pass_av_uvh, g%Domain, clock=id_clock_pass)
817  else
818  call do_group_pass(cs%pass_av_uvh, g%domain, clock=id_clock_pass)
819  endif
820 
821  if (associated(cs%OBC)) then
822  call radiation_open_bdry_conds(cs%OBC, u, u_old_rad_obc, v, v_old_rad_obc, g, us, dt)
823  endif
824 
825 ! h_av = (h_in + h_out)/2 . Going in to this line, h_av = h_in.
826  !$OMP parallel do default(shared)
827  do k=1,nz ; do j=js-2,je+2 ; do i=is-2,ie+2
828  h_av(i,j,k) = 0.5*(h_av(i,j,k) + h(i,j,k))
829  enddo ; enddo ; enddo
830 
831  if (g%nonblocking_updates) &
832  call complete_group_pass(cs%pass_av_uvh, g%Domain, clock=id_clock_pass)
833 
834  !$OMP parallel do default(shared)
835  do k=1,nz
836  do j=js-2,je+2 ; do i=isq-2,ieq+2
837  uhtr(i,j,k) = uhtr(i,j,k) + uh(i,j,k)*dt
838  enddo ; enddo
839  do j=jsq-2,jeq+2 ; do i=is-2,ie+2
840  vhtr(i,j,k) = vhtr(i,j,k) + vh(i,j,k)*dt
841  enddo ; enddo
842  enddo
843 
844  ! The time-averaged free surface height has already been set by the last
845  ! call to btstep.
846 
847  ! Here various terms used in to update the momentum equations are
848  ! offered for time averaging.
849  if (cs%id_PFu > 0) call post_data(cs%id_PFu, cs%PFu, cs%diag)
850  if (cs%id_PFv > 0) call post_data(cs%id_PFv, cs%PFv, cs%diag)
851  if (cs%id_CAu > 0) call post_data(cs%id_CAu, cs%CAu, cs%diag)
852  if (cs%id_CAv > 0) call post_data(cs%id_CAv, cs%CAv, cs%diag)
853 
854  ! Here the thickness fluxes are offered for time averaging.
855  if (cs%id_uh > 0) call post_data(cs%id_uh , uh, cs%diag)
856  if (cs%id_vh > 0) call post_data(cs%id_vh , vh, cs%diag)
857  if (cs%id_uav > 0) call post_data(cs%id_uav, u_av, cs%diag)
858  if (cs%id_vav > 0) call post_data(cs%id_vav, v_av, cs%diag)
859  if (cs%id_u_BT_accel > 0) call post_data(cs%id_u_BT_accel, cs%u_accel_bt, cs%diag)
860  if (cs%id_v_BT_accel > 0) call post_data(cs%id_v_BT_accel, cs%v_accel_bt, cs%diag)
861 
862  if (cs%debug) then
863  call mom_state_chksum("Corrector ", u, v, h, uh, vh, g, gv, us, symmetric=sym)
864  call uvchksum("Corrector avg [uv]", u_av, v_av, g%HI, haloshift=1, symmetric=sym, scale=us%L_T_to_m_s)
865  call hchksum(h_av, "Corrector avg h", g%HI, haloshift=1, scale=gv%H_to_m)
866  ! call MOM_state_chksum("Corrector avg ", u_av, v_av, h_av, uh, vh, G, GV, US)
867  endif
868 
869  if (showcalltree) call calltree_leave("step_MOM_dyn_split_RK2()")
870 

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_test_extern_h(), 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().

Here is the call graph for this function:

Variable Documentation

◆ id_clock_btcalc

integer mom_dynamics_split_rk2::id_clock_btcalc
private

CPU time clock IDs.

Definition at line 228 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

CPU time clock IDs.

Definition at line 228 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

CPU time clock IDs.

Definition at line 228 of file MOM_dynamics_split_RK2.F90.

228 integer :: id_clock_btstep, id_clock_btcalc, id_clock_btforce

Referenced by initialize_dyn_split_rk2(), and step_mom_dyn_split_rk2().

◆ id_clock_continuity

integer mom_dynamics_split_rk2::id_clock_continuity
private

CPU time clock IDs.

Definition at line 227 of file MOM_dynamics_split_RK2.F90.

227 integer :: id_clock_continuity, id_clock_thick_diff

Referenced by initialize_dyn_split_rk2(), and step_mom_dyn_split_rk2().

◆ id_clock_cor

integer mom_dynamics_split_rk2::id_clock_cor

CPU time clock IDs.

Definition at line 225 of file MOM_dynamics_split_RK2.F90.

225 integer :: id_clock_Cor, id_clock_pres, id_clock_vertvisc

Referenced by initialize_dyn_split_rk2(), and step_mom_dyn_split_rk2().

◆ id_clock_horvisc

integer mom_dynamics_split_rk2::id_clock_horvisc
private

CPU time clock IDs.

Definition at line 226 of file MOM_dynamics_split_RK2.F90.

226 integer :: id_clock_horvisc, id_clock_mom_update

Referenced by initialize_dyn_split_rk2(), and step_mom_dyn_split_rk2().

◆ id_clock_mom_update

integer mom_dynamics_split_rk2::id_clock_mom_update
private

CPU time clock IDs.

Definition at line 226 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

CPU time clock IDs.

Definition at line 229 of file MOM_dynamics_split_RK2.F90.

229 integer :: id_clock_pass, id_clock_pass_init

Referenced by initialize_dyn_split_rk2(), and step_mom_dyn_split_rk2().

◆ id_clock_pass_init

integer mom_dynamics_split_rk2::id_clock_pass_init
private

CPU time clock IDs.

Definition at line 229 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

CPU time clock IDs.

Definition at line 225 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

CPU time clock IDs.

Definition at line 227 of file MOM_dynamics_split_RK2.F90.

◆ id_clock_vertvisc

integer mom_dynamics_split_rk2::id_clock_vertvisc
private

CPU time clock IDs.

Definition at line 225 of file MOM_dynamics_split_RK2.F90.

Referenced by initialize_dyn_split_rk2(), and step_mom_dyn_split_rk2().