MOM6
Data Types | Functions/Subroutines
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...
 

Functions/Subroutines

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, cont_stencil)
 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 1443 of file MOM_dynamics_split_RK2.F90.

1444  type(MOM_dyn_split_RK2_CS), pointer :: CS !< module control structure
1445 
1446  dealloc_(cs%diffu) ; dealloc_(cs%diffv)
1447  dealloc_(cs%CAu) ; dealloc_(cs%CAv)
1448  dealloc_(cs%PFu) ; dealloc_(cs%PFv)
1449 
1450  if (associated(cs%taux_bot)) deallocate(cs%taux_bot)
1451  if (associated(cs%tauy_bot)) deallocate(cs%tauy_bot)
1452  dealloc_(cs%uhbt) ; dealloc_(cs%vhbt)
1453  dealloc_(cs%u_accel_bt) ; dealloc_(cs%v_accel_bt)
1454  dealloc_(cs%visc_rem_u) ; dealloc_(cs%visc_rem_v)
1455 
1456  dealloc_(cs%eta) ; dealloc_(cs%eta_PF) ; dealloc_(cs%pbce)
1457  dealloc_(cs%h_av) ; dealloc_(cs%u_av) ; dealloc_(cs%v_av)
1458 
1459  call dealloc_bt_cont_type(cs%BT_cont)
1460 
1461  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,
integer, intent(out), optional  cont_stencil 
)

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
[out]cont_stencilThe stencil for thickness from the continuity solver.

Definition at line 1072 of file MOM_dynamics_split_RK2.F90.

1077  type(ocean_grid_type), intent(inout) :: G !< ocean grid structure
1078  type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
1079  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1080  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
1081  intent(inout) :: u !< zonal velocity [L T-1 ~> m s-1]
1082  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1083  intent(inout) :: v !< merid velocity [L T-1 ~> m s-1]
1084  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) , intent(inout) :: h !< layer thickness [H ~> m or kg m-2]
1085  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
1086  target, intent(inout) :: uh !< zonal volume/mass transport [H L2 T-1 ~> m3 s-1 or kg s-1]
1087  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1088  target, intent(inout) :: vh !< merid volume/mass transport [H L2 T-1 ~> m3 s-1 or kg s-1]
1089  real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: eta !< free surface height or column mass [H ~> m or kg m-2]
1090  type(time_type), target, intent(in) :: Time !< current model time
1091  type(param_file_type), intent(in) :: param_file !< parameter file for parsing
1092  type(diag_ctrl), target, intent(inout) :: diag !< to control diagnostics
1093  type(MOM_dyn_split_RK2_CS), pointer :: CS !< module control structure
1094  type(MOM_restart_CS), pointer :: restart_CS !< restart control structure
1095  real, intent(in) :: dt !< time step [T ~> s]
1096  type(accel_diag_ptrs), target, intent(inout) :: Accel_diag !< points to momentum equation terms for
1097  !! budget analysis
1098  type(cont_diag_ptrs), target, intent(inout) :: Cont_diag !< points to terms in continuity equation
1099  type(ocean_internal_state), intent(inout) :: MIS !< "MOM6 internal state" used to pass
1100  !! diagnostic pointers
1101  type(VarMix_CS), pointer :: VarMix !< points to spatially variable viscosities
1102  type(MEKE_type), pointer :: MEKE !< points to mesoscale eddy kinetic energy fields
1103 ! type(Barotropic_CS), pointer :: Barotropic_CSp !< Pointer to the control structure for
1104 ! !! the barotropic module
1105  type(thickness_diffuse_CS), pointer :: thickness_diffuse_CSp !< Pointer to the control structure
1106  !! used for the isopycnal height diffusive transport.
1107  type(ocean_OBC_type), pointer :: OBC !< points to OBC related fields
1108  type(update_OBC_CS), pointer :: update_OBC_CSp !< points to OBC update related fields
1109  type(ALE_CS), pointer :: ALE_CSp !< points to ALE control structure
1110  type(set_visc_CS), pointer :: setVisc_CSp !< points to the set_visc control structure.
1111  type(vertvisc_type), intent(inout) :: visc !< vertical viscosities, bottom drag, and related
1112  type(directories), intent(in) :: dirs !< contains directory paths
1113  integer, target, intent(inout) :: ntrunc !< A target for the variable that records
1114  !! the number of times the velocity is
1115  !! truncated (this should be 0).
1116  logical, intent(out) :: calc_dtbt !< If true, recalculate the barotropic time step
1117  integer, optional, intent(out) :: cont_stencil !< The stencil for thickness
1118  !! from the continuity solver.
1119 
1120  ! local variables
1121  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_tmp
1122  character(len=40) :: mdl = "MOM_dynamics_split_RK2" ! This module's name.
1123  ! This include declares and sets the variable "version".
1124 # include "version_variable.h"
1125  character(len=48) :: thickness_units, flux_units, eta_rest_name
1126  real :: H_rescale ! A rescaling factor for thicknesses from the representation in
1127  ! a restart file to the internal representation in this run.
1128  real :: vel_rescale ! A rescaling factor for velocities from the representation in
1129  ! a restart file to the internal representation in this run.
1130  real :: uH_rescale ! A rescaling factor for thickness transports from the representation in
1131  ! a restart file to the internal representation in this run.
1132  real :: accel_rescale ! A rescaling factor for accelerations from the representation in
1133  ! a restart file to the internal representation in this run.
1134  real :: H_convert
1135  type(group_pass_type) :: pass_av_h_uvh
1136  logical :: use_tides, debug_truncations
1137 
1138  integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz
1139  integer :: IsdB, IedB, JsdB, JedB
1140  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
1141  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
1142  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
1143 
1144  if (.not.associated(cs)) call mom_error(fatal, &
1145  "initialize_dyn_split_RK2 called with an unassociated control structure.")
1146  if (cs%module_is_initialized) then
1147  call mom_error(warning, "initialize_dyn_split_RK2 called with a control "// &
1148  "structure that has already been initialized.")
1149  return
1150  endif
1151  cs%module_is_initialized = .true.
1152 
1153  cs%diag => diag
1154 
1155  call log_version(param_file, mdl, version, "")
1156  call get_param(param_file, mdl, "TIDES", use_tides, &
1157  "If true, apply tidal momentum forcing.", default=.false.)
1158  call get_param(param_file, mdl, "BE", cs%be, &
1159  "If SPLIT is true, BE determines the relative weighting "//&
1160  "of a 2nd-order Runga-Kutta baroclinic time stepping "//&
1161  "scheme (0.5) and a backward Euler scheme (1) that is "//&
1162  "used for the Coriolis and inertial terms. BE may be "//&
1163  "from 0.5 to 1, but instability may occur near 0.5. "//&
1164  "BE is also applicable if SPLIT is false and USE_RK2 "//&
1165  "is true.", units="nondim", default=0.6)
1166  call get_param(param_file, mdl, "BEGW", cs%begw, &
1167  "If SPLIT is true, BEGW is a number from 0 to 1 that "//&
1168  "controls the extent to which the treatment of gravity "//&
1169  "waves is forward-backward (0) or simulated backward "//&
1170  "Euler (1). 0 is almost always used. "//&
1171  "If SPLIT is false and USE_RK2 is true, BEGW can be "//&
1172  "between 0 and 0.5 to damp gravity waves.", &
1173  units="nondim", default=0.0)
1174 
1175  call get_param(param_file, mdl, "SPLIT_BOTTOM_STRESS", cs%split_bottom_stress, &
1176  "If true, provide the bottom stress calculated by the "//&
1177  "vertical viscosity to the barotropic solver.", default=.false.)
1178  call get_param(param_file, mdl, "BT_USE_LAYER_FLUXES", cs%BT_use_layer_fluxes, &
1179  "If true, use the summed layered fluxes plus an "//&
1180  "adjustment due to the change in the barotropic velocity "//&
1181  "in the barotropic continuity equation.", default=.true.)
1182  call get_param(param_file, mdl, "DEBUG", cs%debug, &
1183  "If true, write out verbose debugging data.", &
1184  default=.false., debuggingparam=.true.)
1185  call get_param(param_file, mdl, "DEBUG_OBC", cs%debug_OBC, default=.false.)
1186  call get_param(param_file, mdl, "DEBUG_TRUNCATIONS", debug_truncations, &
1187  default=.false.)
1188 
1189  allocate(cs%taux_bot(isdb:iedb,jsd:jed)) ; cs%taux_bot(:,:) = 0.0
1190  allocate(cs%tauy_bot(isd:ied,jsdb:jedb)) ; cs%tauy_bot(:,:) = 0.0
1191 
1192  alloc_(cs%uhbt(isdb:iedb,jsd:jed)) ; cs%uhbt(:,:) = 0.0
1193  alloc_(cs%vhbt(isd:ied,jsdb:jedb)) ; cs%vhbt(:,:) = 0.0
1194  alloc_(cs%visc_rem_u(isdb:iedb,jsd:jed,nz)) ; cs%visc_rem_u(:,:,:) = 0.0
1195  alloc_(cs%visc_rem_v(isd:ied,jsdb:jedb,nz)) ; cs%visc_rem_v(:,:,:) = 0.0
1196  alloc_(cs%eta_PF(isd:ied,jsd:jed)) ; cs%eta_PF(:,:) = 0.0
1197  alloc_(cs%pbce(isd:ied,jsd:jed,nz)) ; cs%pbce(:,:,:) = 0.0
1198 
1199  alloc_(cs%u_accel_bt(isdb:iedb,jsd:jed,nz)) ; cs%u_accel_bt(:,:,:) = 0.0
1200  alloc_(cs%v_accel_bt(isd:ied,jsdb:jedb,nz)) ; cs%v_accel_bt(:,:,:) = 0.0
1201 
1202  mis%diffu => cs%diffu
1203  mis%diffv => cs%diffv
1204  mis%PFu => cs%PFu
1205  mis%PFv => cs%PFv
1206  mis%CAu => cs%CAu
1207  mis%CAv => cs%CAv
1208  mis%pbce => cs%pbce
1209  mis%u_accel_bt => cs%u_accel_bt
1210  mis%v_accel_bt => cs%v_accel_bt
1211  mis%u_av => cs%u_av
1212  mis%v_av => cs%v_av
1213 
1214  cs%ADp => accel_diag
1215  cs%CDp => cont_diag
1216  accel_diag%diffu => cs%diffu
1217  accel_diag%diffv => cs%diffv
1218  accel_diag%PFu => cs%PFu
1219  accel_diag%PFv => cs%PFv
1220  accel_diag%CAu => cs%CAu
1221  accel_diag%CAv => cs%CAv
1222 
1223 ! Accel_diag%pbce => CS%pbce
1224 ! Accel_diag%u_accel_bt => CS%u_accel_bt ; Accel_diag%v_accel_bt => CS%v_accel_bt
1225 ! Accel_diag%u_av => CS%u_av ; Accel_diag%v_av => CS%v_av
1226 
1227  id_clock_pass_init = cpu_clock_id('(Ocean init message passing)', &
1228  grain=clock_routine)
1229 
1230  call continuity_init(time, g, gv, us, param_file, diag, cs%continuity_CSp)
1231  if (present(cont_stencil)) cont_stencil = continuity_stencil(cs%continuity_CSp)
1232  call coriolisadv_init(time, g, gv, us, param_file, diag, cs%ADp, cs%CoriolisAdv_CSp)
1233  if (use_tides) call tidal_forcing_init(time, g, param_file, cs%tides_CSp)
1234  call pressureforce_init(time, g, gv, us, param_file, diag, cs%PressureForce_CSp, &
1235  cs%tides_CSp)
1236  call hor_visc_init(time, g, us, param_file, diag, cs%hor_visc_CSp, meke, adp=cs%ADp)
1237  call vertvisc_init(mis, time, g, gv, us, param_file, diag, cs%ADp, dirs, &
1238  ntrunc, cs%vertvisc_CSp)
1239  if (.not.associated(setvisc_csp)) call mom_error(fatal, &
1240  "initialize_dyn_split_RK2 called with setVisc_CSp unassociated.")
1241  cs%set_visc_CSp => setvisc_csp
1242  call updatecfltruncationvalue(time, cs%vertvisc_CSp, &
1243  activate=is_new_run(restart_cs) )
1244 
1245  if (associated(ale_csp)) cs%ALE_CSp => ale_csp
1246  if (associated(obc)) then
1247  cs%OBC => obc
1248  if (obc%ramp) call update_obc_ramp(time, cs%OBC, &
1249  activate=is_new_run(restart_cs) )
1250  endif
1251  if (associated(update_obc_csp)) cs%update_OBC_CSp => update_obc_csp
1252 
1253  eta_rest_name = "sfc" ; if (.not.gv%Boussinesq) eta_rest_name = "p_bot"
1254  if (.not. query_initialized(cs%eta, trim(eta_rest_name), restart_cs)) then
1255  ! Estimate eta based on the layer thicknesses - h. With the Boussinesq
1256  ! approximation, eta is the free surface height anomaly, while without it
1257  ! eta is the mass of ocean per unit area. eta always has the same
1258  ! dimensions as h, either m or kg m-3.
1259  ! CS%eta(:,:) = 0.0 already from initialization.
1260  if (gv%Boussinesq) then
1261  do j=js,je ; do i=is,ie ; cs%eta(i,j) = -gv%Z_to_H * g%bathyT(i,j) ; enddo ; enddo
1262  endif
1263  do k=1,nz ; do j=js,je ; do i=is,ie
1264  cs%eta(i,j) = cs%eta(i,j) + h(i,j,k)
1265  enddo ; enddo ; enddo
1266  elseif ((gv%m_to_H_restart /= 0.0) .and. (gv%m_to_H_restart /= gv%m_to_H)) then
1267  h_rescale = gv%m_to_H / gv%m_to_H_restart
1268  do j=js,je ; do i=is,ie ; cs%eta(i,j) = h_rescale * cs%eta(i,j) ; enddo ; enddo
1269  endif
1270  ! Copy eta into an output array.
1271  do j=js,je ; do i=is,ie ; eta(i,j) = cs%eta(i,j) ; enddo ; enddo
1272 
1273  call barotropic_init(u, v, h, cs%eta, time, g, gv, us, param_file, diag, &
1274  cs%barotropic_CSp, restart_cs, calc_dtbt, cs%BT_cont, &
1275  cs%tides_CSp)
1276 
1277  if (.not. query_initialized(cs%diffu,"diffu",restart_cs) .or. &
1278  .not. query_initialized(cs%diffv,"diffv",restart_cs)) then
1279  call horizontal_viscosity(u, v, h, cs%diffu, cs%diffv, meke, varmix, &
1280  g, gv, us, cs%hor_visc_CSp, &
1281  obc=cs%OBC, bt=cs%barotropic_CSp, &
1282  td=thickness_diffuse_csp)
1283  else
1284  if ( (us%s_to_T_restart * us%m_to_L_restart /= 0.0) .and. &
1285  (us%m_to_L * us%s_to_T_restart**2 /= us%m_to_L_restart * us%s_to_T**2) ) then
1286  accel_rescale = (us%m_to_L * us%s_to_T_restart**2) / (us%m_to_L_restart * us%s_to_T**2)
1287  do k=1,nz ; do j=js,je ; do i=g%IscB,g%IecB
1288  cs%diffu(i,j,k) = accel_rescale * cs%diffu(i,j,k)
1289  enddo ; enddo ; enddo
1290  do k=1,nz ; do j=g%JscB,g%JecB ; do i=is,ie
1291  cs%diffv(i,j,k) = accel_rescale * cs%diffv(i,j,k)
1292  enddo ; enddo ; enddo
1293  endif
1294  endif
1295 
1296  if (.not. query_initialized(cs%u_av,"u2", restart_cs) .or. &
1297  .not. query_initialized(cs%u_av,"v2", restart_cs)) then
1298  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
1299  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
1300  elseif ( (us%s_to_T_restart * us%m_to_L_restart /= 0.0) .and. &
1301  ((us%m_to_L * us%s_to_T_restart) /= (us%m_to_L_restart * us%s_to_T)) ) then
1302  vel_rescale = (us%m_to_L * us%s_to_T_restart) / (us%m_to_L_restart * us%s_to_T)
1303  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
1304  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
1305  endif
1306 
1307  ! This call is just here to initialize uh and vh.
1308  if (.not. query_initialized(uh,"uh",restart_cs) .or. &
1309  .not. query_initialized(vh,"vh",restart_cs)) then
1310  do k=1,nz ; do j=jsd,jed ; do i=isd,ied ; h_tmp(i,j,k) = h(i,j,k) ; enddo ; enddo ; enddo
1311  call continuity(u, v, h, h_tmp, uh, vh, dt, g, gv, us, cs%continuity_CSp, obc=cs%OBC)
1312  call pass_var(h_tmp, g%Domain, clock=id_clock_pass_init)
1313  do k=1,nz ; do j=jsd,jed ; do i=isd,ied
1314  cs%h_av(i,j,k) = 0.5*(h(i,j,k) + h_tmp(i,j,k))
1315  enddo ; enddo ; enddo
1316  else
1317  if (.not. query_initialized(cs%h_av,"h2",restart_cs)) then
1318  cs%h_av(:,:,:) = h(:,:,:)
1319  elseif ((gv%m_to_H_restart /= 0.0) .and. (gv%m_to_H_restart /= gv%m_to_H)) then
1320  h_rescale = gv%m_to_H / gv%m_to_H_restart
1321  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
1322  endif
1323  if ( (gv%m_to_H_restart * us%s_to_T_restart * us%m_to_L_restart /= 0.0) .and. &
1324  ((gv%m_to_H * us%m_to_L**2 * us%s_to_T_restart) /= &
1325  (gv%m_to_H_restart * us%m_to_L_restart**2 * us%s_to_T)) ) then
1326  uh_rescale = (gv%m_to_H * us%m_to_L**2 * us%s_to_T_restart) / &
1327  (gv%m_to_H_restart * us%m_to_L_restart**2 * us%s_to_T)
1328  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
1329  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
1330  endif
1331  endif
1332 
1333  call cpu_clock_begin(id_clock_pass_init)
1334  call create_group_pass(pass_av_h_uvh, cs%u_av, cs%v_av, g%Domain, halo=2)
1335  call create_group_pass(pass_av_h_uvh, cs%h_av, g%Domain, halo=2)
1336  call create_group_pass(pass_av_h_uvh, uh, vh, g%Domain, halo=2)
1337  call do_group_pass(pass_av_h_uvh, g%Domain)
1338  call cpu_clock_end(id_clock_pass_init)
1339 
1340  flux_units = get_flux_units(gv)
1341  h_convert = gv%H_to_m ; if (.not.gv%Boussinesq) h_convert = gv%H_to_kg_m2
1342  cs%id_uh = register_diag_field('ocean_model', 'uh', diag%axesCuL, time, &
1343  'Zonal Thickness Flux', flux_units, y_cell_method='sum', v_extensive=.true., &
1344  conversion=h_convert*us%L_to_m**2*us%s_to_T)
1345  cs%id_vh = register_diag_field('ocean_model', 'vh', diag%axesCvL, time, &
1346  'Meridional Thickness Flux', flux_units, x_cell_method='sum', v_extensive=.true., &
1347  conversion=h_convert*us%L_to_m**2*us%s_to_T)
1348 
1349  cs%id_CAu = register_diag_field('ocean_model', 'CAu', diag%axesCuL, time, &
1350  'Zonal Coriolis and Advective Acceleration', 'm s-2', conversion=us%L_T2_to_m_s2)
1351  cs%id_CAv = register_diag_field('ocean_model', 'CAv', diag%axesCvL, time, &
1352  'Meridional Coriolis and Advective Acceleration', 'm s-2', conversion=us%L_T2_to_m_s2)
1353  cs%id_PFu = register_diag_field('ocean_model', 'PFu', diag%axesCuL, time, &
1354  'Zonal Pressure Force Acceleration', 'm s-2', conversion=us%L_T2_to_m_s2)
1355  cs%id_PFv = register_diag_field('ocean_model', 'PFv', diag%axesCvL, time, &
1356  'Meridional Pressure Force Acceleration', 'm s-2', conversion=us%L_T2_to_m_s2)
1357 
1358  !CS%id_hf_PFu = register_diag_field('ocean_model', 'hf_PFu', diag%axesCuL, Time, &
1359  ! 'Fractional Thickness-weighted Zonal Pressure Force Acceleration', 'm s-2', &
1360  ! v_extensive=.true., conversion=US%L_T2_to_m_s2)
1361  !if(CS%id_hf_PFu > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_u,IsdB,IedB,jsd,jed,nz)
1362 
1363  !CS%id_hf_PFv = register_diag_field('ocean_model', 'hf_PFv', diag%axesCvL, Time, &
1364  ! 'Fractional Thickness-weighted Meridional Pressure Force Acceleration', 'm s-2', &
1365  ! v_extensive=.true., conversion=US%L_T2_to_m_s2)
1366  !if(CS%id_hf_PFv > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
1367 
1368  !CS%id_hf_CAu = register_diag_field('ocean_model', 'hf_CAu', diag%axesCuL, Time, &
1369  ! 'Fractional Thickness-weighted Zonal Coriolis and Advective Acceleration', 'm s-2', &
1370  ! v_extensive=.true., conversion=US%L_T2_to_m_s2)
1371  !if(CS%id_hf_CAu > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_u,IsdB,IedB,jsd,jed,nz)
1372 
1373  !CS%id_hf_CAv = register_diag_field('ocean_model', 'hf_CAv', diag%axesCvL, Time, &
1374  ! 'Fractional Thickness-weighted Meridional Coriolis and Advective Acceleration', 'm s-2', &
1375  ! v_extensive=.true., conversion=US%L_T2_to_m_s2)
1376  !if(CS%id_hf_CAv > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
1377 
1378  cs%id_hf_PFu_2d = register_diag_field('ocean_model', 'hf_PFu_2d', diag%axesCu1, time, &
1379  'Depth-sum Fractional Thickness-weighted Zonal Pressure Force Acceleration', 'm s-2', &
1380  conversion=us%L_T2_to_m_s2)
1381  if(cs%id_hf_PFu_2d > 0) call safe_alloc_ptr(cs%ADp%diag_hfrac_u,isdb,iedb,jsd,jed,nz)
1382 
1383  cs%id_hf_PFv_2d = register_diag_field('ocean_model', 'hf_PFv_2d', diag%axesCv1, time, &
1384  'Depth-sum Fractional Thickness-weighted Meridional Pressure Force Acceleration', 'm s-2', &
1385  conversion=us%L_T2_to_m_s2)
1386  if(cs%id_hf_PFv_2d > 0) call safe_alloc_ptr(cs%ADp%diag_hfrac_v,isd,ied,jsd,jedb,nz)
1387 
1388  cs%id_hf_CAu_2d = register_diag_field('ocean_model', 'hf_CAu_2d', diag%axesCu1, time, &
1389  'Depth-sum Fractional Thickness-weighted Zonal Coriolis and Advective Acceleration', 'm s-2', &
1390  conversion=us%L_T2_to_m_s2)
1391  if(cs%id_hf_CAu_2d > 0) call safe_alloc_ptr(cs%ADp%diag_hfrac_u,isdb,iedb,jsd,jed,nz)
1392 
1393  cs%id_hf_CAv_2d = register_diag_field('ocean_model', 'hf_CAv_2d', diag%axesCv1, time, &
1394  'Depth-sum Fractional Thickness-weighted Meridional Coriolis and Advective Acceleration', 'm s-2', &
1395  conversion=us%L_T2_to_m_s2)
1396  if(cs%id_hf_CAv_2d > 0) call safe_alloc_ptr(cs%ADp%diag_hfrac_v,isd,ied,jsd,jedb,nz)
1397 
1398  cs%id_uav = register_diag_field('ocean_model', 'uav', diag%axesCuL, time, &
1399  'Barotropic-step Averaged Zonal Velocity', 'm s-1', conversion=us%L_T_to_m_s)
1400  cs%id_vav = register_diag_field('ocean_model', 'vav', diag%axesCvL, time, &
1401  'Barotropic-step Averaged Meridional Velocity', 'm s-1', conversion=us%L_T_to_m_s)
1402 
1403  cs%id_u_BT_accel = register_diag_field('ocean_model', 'u_BT_accel', diag%axesCuL, time, &
1404  'Barotropic Anomaly Zonal Acceleration', 'm s-2', conversion=us%L_T2_to_m_s2)
1405  cs%id_v_BT_accel = register_diag_field('ocean_model', 'v_BT_accel', diag%axesCvL, time, &
1406  'Barotropic Anomaly Meridional Acceleration', 'm s-2', conversion=us%L_T2_to_m_s2)
1407 
1408  !CS%id_hf_u_BT_accel = register_diag_field('ocean_model', 'hf_u_BT_accel', diag%axesCuL, Time, &
1409  ! 'Fractional Thickness-weighted Barotropic Anomaly Zonal Acceleration', 'm s-2', &
1410  ! v_extensive=.true., conversion=US%L_T2_to_m_s2)
1411  !if(CS%id_hf_u_BT_accel > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_u,IsdB,IedB,jsd,jed,nz)
1412 
1413  !CS%id_hf_v_BT_accel = register_diag_field('ocean_model', 'hf_v_BT_accel', diag%axesCvL, Time, &
1414  ! 'Fractional Thickness-weighted Barotropic Anomaly Meridional Acceleration', 'm s-2', &
1415  ! v_extensive=.true., conversion=US%L_T2_to_m_s2)
1416  !if(CS%id_hf_v_BT_accel > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
1417 
1418  cs%id_hf_u_BT_accel_2d = register_diag_field('ocean_model', 'hf_u_BT_accel_2d', diag%axesCu1, time, &
1419  'Depth-sum Fractional Thickness-weighted Barotropic Anomaly Zonal Acceleration', 'm s-2', &
1420  conversion=us%L_T2_to_m_s2)
1421  if(cs%id_hf_u_BT_accel_2d > 0) call safe_alloc_ptr(cs%ADp%diag_hfrac_u,isdb,iedb,jsd,jed,nz)
1422 
1423  cs%id_hf_v_BT_accel_2d = register_diag_field('ocean_model', 'hf_v_BT_accel_2d', diag%axesCv1, time, &
1424  'Depth-sum Fractional Thickness-weighted Barotropic Anomaly Meridional Acceleration', 'm s-2', &
1425  conversion=us%L_T2_to_m_s2)
1426  if(cs%id_hf_v_BT_accel_2d > 0) call safe_alloc_ptr(cs%ADp%diag_hfrac_v,isd,ied,jsd,jedb,nz)
1427 
1428  id_clock_cor = cpu_clock_id('(Ocean Coriolis & mom advection)', grain=clock_module)
1429  id_clock_continuity = cpu_clock_id('(Ocean continuity equation)', grain=clock_module)
1430  id_clock_pres = cpu_clock_id('(Ocean pressure force)', grain=clock_module)
1431  id_clock_vertvisc = cpu_clock_id('(Ocean vertical viscosity)', grain=clock_module)
1432  id_clock_horvisc = cpu_clock_id('(Ocean horizontal viscosity)', grain=clock_module)
1433  id_clock_mom_update = cpu_clock_id('(Ocean momentum increments)', grain=clock_module)
1434  id_clock_pass = cpu_clock_id('(Ocean message passing)', grain=clock_module)
1435  id_clock_btcalc = cpu_clock_id('(Ocean barotropic mode calc)', grain=clock_module)
1436  id_clock_btstep = cpu_clock_id('(Ocean barotropic mode stepping)', grain=clock_module)
1437  id_clock_btforce = cpu_clock_id('(Ocean barotropic forcing calc)', grain=clock_module)
1438 

◆ 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 1000 of file MOM_dynamics_split_RK2.F90.

1001  type(hor_index_type), intent(in) :: HI !< Horizontal index structure
1002  type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
1003  type(param_file_type), intent(in) :: param_file !< parameter file
1004  type(MOM_dyn_split_RK2_CS), pointer :: CS !< module control structure
1005  type(MOM_restart_CS), pointer :: restart_CS !< restart control structure
1006  real, dimension(SZIB_(HI),SZJ_(HI),SZK_(GV)), &
1007  target, intent(inout) :: uh !< zonal volume/mass transport [H L2 T-1 ~> m3 s-1 or kg s-1]
1008  real, dimension(SZI_(HI),SZJB_(HI),SZK_(GV)), &
1009  target, intent(inout) :: vh !< merid volume/mass transport [H L2 T-1 ~> m3 s-1 or kg s-1]
1010 
1011  type(vardesc) :: vd(2)
1012  character(len=40) :: mdl = "MOM_dynamics_split_RK2" ! This module's name.
1013  character(len=48) :: thickness_units, flux_units
1014 
1015  integer :: isd, ied, jsd, jed, nz, IsdB, IedB, JsdB, JedB
1016 
1017  isd = hi%isd ; ied = hi%ied ; jsd = hi%jsd ; jed = hi%jed ; nz = gv%ke
1018  isdb = hi%IsdB ; iedb = hi%IedB ; jsdb = hi%JsdB ; jedb = hi%JedB
1019 
1020  ! This is where a control structure specific to this module would be allocated.
1021  if (associated(cs)) then
1022  call mom_error(warning, "register_restarts_dyn_split_RK2 called with an associated "// &
1023  "control structure.")
1024  return
1025  endif
1026  allocate(cs)
1027 
1028  alloc_(cs%diffu(isdb:iedb,jsd:jed,nz)) ; cs%diffu(:,:,:) = 0.0
1029  alloc_(cs%diffv(isd:ied,jsdb:jedb,nz)) ; cs%diffv(:,:,:) = 0.0
1030  alloc_(cs%CAu(isdb:iedb,jsd:jed,nz)) ; cs%CAu(:,:,:) = 0.0
1031  alloc_(cs%CAv(isd:ied,jsdb:jedb,nz)) ; cs%CAv(:,:,:) = 0.0
1032  alloc_(cs%PFu(isdb:iedb,jsd:jed,nz)) ; cs%PFu(:,:,:) = 0.0
1033  alloc_(cs%PFv(isd:ied,jsdb:jedb,nz)) ; cs%PFv(:,:,:) = 0.0
1034 
1035  alloc_(cs%eta(isd:ied,jsd:jed)) ; cs%eta(:,:) = 0.0
1036  alloc_(cs%u_av(isdb:iedb,jsd:jed,nz)) ; cs%u_av(:,:,:) = 0.0
1037  alloc_(cs%v_av(isd:ied,jsdb:jedb,nz)) ; cs%v_av(:,:,:) = 0.0
1038  alloc_(cs%h_av(isd:ied,jsd:jed,nz)) ; cs%h_av(:,:,:) = gv%Angstrom_H
1039 
1040  thickness_units = get_thickness_units(gv)
1041  flux_units = get_flux_units(gv)
1042 
1043  if (gv%Boussinesq) then
1044  vd(1) = var_desc("sfc",thickness_units,"Free surface Height",'h','1')
1045  else
1046  vd(1) = var_desc("p_bot",thickness_units,"Bottom Pressure",'h','1')
1047  endif
1048  call register_restart_field(cs%eta, vd(1), .false., restart_cs)
1049 
1050  vd(1) = var_desc("u2","m s-1","Auxiliary Zonal velocity",'u','L')
1051  vd(2) = var_desc("v2","m s-1","Auxiliary Meridional velocity",'v','L')
1052  call register_restart_pair(cs%u_av, cs%v_av, vd(1), vd(2), .false., restart_cs)
1053 
1054  vd(1) = var_desc("h2",thickness_units,"Auxiliary Layer Thickness",'h','L')
1055  call register_restart_field(cs%h_av, vd(1), .false., restart_cs)
1056 
1057  vd(1) = var_desc("uh",flux_units,"Zonal thickness flux",'u','L')
1058  vd(2) = var_desc("vh",flux_units,"Meridional thickness flux",'v','L')
1059  call register_restart_pair(uh, vh, vd(1), vd(2), .false., restart_cs)
1060 
1061  vd(1) = var_desc("diffu","m s-2","Zonal horizontal viscous acceleration",'u','L')
1062  vd(2) = var_desc("diffv","m s-2","Meridional horizontal viscous acceleration",'v','L')
1063  call register_restart_pair(cs%diffu, cs%diffv, vd(1), vd(2), .false., restart_cs)
1064 
1065  call register_barotropic_restarts(hi, gv, param_file, cs%barotropic_CSp, &
1066  restart_cs)
1067 

◆ 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 the start of this dynamic time step [R L2 T-2 ~> Pa]
p_surf_endsurf pressure at the end of this dynamic time step [R L2 T-2 ~> 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 242 of file MOM_dynamics_split_RK2.F90.

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