MOM_dynamics_split_RK2b.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.

2! See the LICENSE file for licensing information.

3! SPDX-License-Identifier: Apache-2.0

4

5!> Time step the adiabatic dynamic core of MOM using RK2 method with greater use of the

6!! time-filtered velocities and less inheritance of tedencies from the previous step in the

7!! predictor step than in the original MOM_dyanmics_split_RK2.

8module mom_dynamics_split_rk2b

9

10use mom_variables, only : vertvisc_type, thermo_var_ptrs, porous_barrier_type

11use mom_variables, only : bt_cont_type, alloc_bt_cont_type, dealloc_bt_cont_type

12use mom_variables, only : accel_diag_ptrs, ocean_internal_state, cont_diag_ptrs

13use mom_forcing_type, only : mech_forcing

14

15use mom_checksum_packages, only : mom_thermo_chksum, mom_state_chksum, mom_accel_chksum

16use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end

17use mom_cpu_clock, only : clock_component, clock_subcomponent

18use mom_cpu_clock, only : clock_module_driver, clock_module, clock_routine

19use mom_diag_mediator, only : diag_mediator_init, enable_averages

20use mom_diag_mediator, only : disable_averaging, post_data, safe_alloc_ptr

21use mom_diag_mediator, only : post_product_u, post_product_sum_u

22use mom_diag_mediator, only : post_product_v, post_product_sum_v

23use mom_diag_mediator, only : register_diag_field, register_static_field

24use mom_diag_mediator, only : set_diag_mediator_grid, diag_ctrl, diag_update_remap_grids

25use mom_domains, only : to_south, to_west, to_all, cgrid_ne, scalar_pair

26use mom_domains, only : to_north, to_east, omit_corners

27use mom_domains, only : create_group_pass, do_group_pass, group_pass_type

28use mom_domains, only : start_group_pass, complete_group_pass, pass_var, pass_vector

29use mom_debugging, only : hchksum, uvchksum, query_debugging_checks

30use mom_error_handler, only : mom_error, mom_mesg, fatal, warning, is_root_pe

31use mom_error_handler, only : mom_set_verbosity, calltree_showquery

32use mom_error_handler, only : calltree_enter, calltree_leave, calltree_waypoint

33use mom_file_parser, only : get_param, log_version, param_file_type

34use mom_get_input, only : directories

35use mom_io, only : vardesc, var_desc, east_face, north_face

36use mom_restart, only : register_restart_field, register_restart_pair

37use mom_restart, only : query_initialized, set_initialized, save_restart

38use mom_restart, only : only_read_from_restarts

39use mom_restart, only : restart_init, is_new_run, mom_restart_cs

40use mom_time_manager, only : time_type, operator(+)

41use mom_time_manager, only : operator(-), operator(>), operator(*), operator(/)

42

43use mom_ale, only : ale_cs, ale_remap_velocities

44use mom_barotropic, only : barotropic_init, btstep, btcalc, bt_mass_source

45use mom_barotropic, only : register_barotropic_restarts, set_dtbt, barotropic_cs

46use mom_barotropic, only : barotropic_end

47use mom_boundary_update, only : update_obc_data, update_obc_cs

48use mom_continuity, only : continuity, continuity_cs

49use mom_continuity, only : continuity_init, continuity_stencil

50use mom_coriolisadv, only : coradcalc, coriolisadv_cs

51use mom_coriolisadv, only : coriolisadv_init, coriolisadv_end, coriolisadv_stencil

52use mom_debugging, only : check_redundant

53use mom_grid, only : ocean_grid_type

54use mom_harmonic_analysis, only : ha_init, harmonic_analysis_cs

55use mom_hor_index, only : hor_index_type

56use mom_hor_visc, only : horizontal_viscosity, hor_visc_cs

57use mom_hor_visc, only : hor_visc_init, hor_visc_end

58use mom_interface_heights, only : thickness_to_dz, find_col_avg_spv

59use mom_lateral_mixing_coeffs, only : varmix_cs

60use mom_meke_types, only : meke_type

61use mom_open_boundary, only : ocean_obc_type, radiation_open_bdry_conds

62use mom_open_boundary, only : open_boundary_zero_normal_flow, open_boundary_query

63use mom_open_boundary, only : open_boundary_test_extern_h, update_obc_ramp

64use mom_open_boundary, only : copy_thickness_reservoirs

65use mom_open_boundary, only : update_segment_thickness_reservoirs

66use mom_pressureforce, only : pressureforce, pressureforce_cs

67use mom_pressureforce, only : pressureforce_init

68use mom_set_visc, only : set_viscous_ml, set_visc_cs

69use mom_thickness_diffuse, only : thickness_diffuse_cs

70use mom_self_attr_load, only : sal_cs

71use mom_self_attr_load, only : sal_init, sal_end

72use mom_tidal_forcing, only : tidal_forcing_cs

73use mom_tidal_forcing, only : tidal_forcing_init, tidal_forcing_end

74use mom_unit_scaling, only : unit_scale_type

75use mom_vert_friction, only : vertvisc, vertvisc_coef, vertvisc_remnant

76use mom_vert_friction, only : vertvisc_init, vertvisc_end, vertvisc_cs

77use mom_vert_friction, only : updatecfltruncationvalue, vertfpmix

78use mom_verticalgrid, only : verticalgrid_type, get_thickness_units

79use mom_verticalgrid, only : get_flux_units, get_tr_flux_units

80use mom_wave_interface, only : wave_parameters_cs, stokes_pgf

81

82implicit none ; private

83

84#include <MOM_memory.h>

85

86!> MOM_dynamics_split_RK2b module control structure

87type, public :: mom_dyn_split_rk2b_cs ; private

88 real allocable_, dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: &

89 cau, & !< CAu = f*v - u.grad(u) [L T-2 ~> m s-2]

90 cau_pred, & !< The predictor step value of CAu = f*v - u.grad(u) [L T-2 ~> m s-2]

91 pfu, & !< PFu = -dM/dx [L T-2 ~> m s-2]

92 pfu_stokes, & !< PFu_Stokes = -d/dx int_r (u_L*duS/dr) [L T-2 ~> m s-2]

93 diffu !< Zonal acceleration due to convergence of the along-isopycnal stress tensor [L T-2 ~> m s-2]

94

95 real allocable_, dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: &

96 cav, & !< CAv = -f*u - u.grad(v) [L T-2 ~> m s-2]

97 cav_pred, & !< The predictor step value of CAv = -f*u - u.grad(v) [L T-2 ~> m s-2]

98 pfv, & !< PFv = -dM/dy [L T-2 ~> m s-2]

99 pfv_stokes, & !< PFv_Stokes = -d/dy int_r (v_L*dvS/dr) [L T-2 ~> m s-2]

100 diffv !< Meridional acceleration due to convergence of the along-isopycnal stress tensor [L T-2 ~> m s-2]

101

102 real allocable_, dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: visc_rem_u

103 !< Both the fraction of the zonal momentum originally in a

104 !! layer that remains after a time-step of viscosity, and the

105 !! fraction of a time-step worth of a barotropic acceleration

106 !! that a layer experiences after viscosity is applied [nondim].

107 !! Nondimensional between 0 (at the bottom) and 1 (far above).

108 real allocable_, dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: u_accel_bt

109 !< The zonal layer accelerations due to the difference between

110 !! the barotropic accelerations and the baroclinic accelerations

111 !! that were fed into the barotopic calculation [L T-2 ~> m s-2]

112 real allocable_, dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: visc_rem_v

113 !< Both the fraction of the meridional momentum originally in

114 !! a layer that remains after a time-step of viscosity, and the

115 !! fraction of a time-step worth of a barotropic acceleration

116 !! that a layer experiences after viscosity is applied [nondim].

117 !! Nondimensional between 0 (at the bottom) and 1 (far above).

118 real allocable_, dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: v_accel_bt

119 !< The meridional layer accelerations due to the difference between

120 !! the barotropic accelerations and the baroclinic accelerations

121 !! that were fed into the barotopic calculation [L T-2 ~> m s-2]

122

123 ! The following variables are only used with the split time stepping scheme.

124 real allocable_, dimension(NIMEM_,NJMEM_) :: eta !< Instantaneous free surface height (in Boussinesq

125 !! mode) or column mass anomaly (in non-Boussinesq

126 !! mode) [H ~> m or kg m-2]

127 real allocable_, dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: u_av !< layer x-velocity with vertical mean replaced by

128 !! time-mean barotropic velocity over a baroclinic

129 !! timestep [L T-1 ~> m s-1]

130 real allocable_, dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: v_av !< layer y-velocity with vertical mean replaced by

131 !! time-mean barotropic velocity over a baroclinic

132 !! timestep [L T-1 ~> m s-1]

133 real allocable_, dimension(NIMEM_,NJMEM_,NKMEM_) :: h_av !< arithmetic mean of two successive layer

134 !! thicknesses [H ~> m or kg m-2]

135 real allocable_, dimension(NIMEM_,NJMEM_) :: eta_pf !< instantaneous SSH used in calculating PFu and

136 !! PFv [H ~> m or kg m-2]

137 real allocable_, dimension(NIMEMB_PTR_,NJMEM_) :: uhbt !< average x-volume or mass flux determined by the

138 !! barotropic solver [H L2 T-1 ~> m3 s-1 or kg s-1].

139 !! uhbt is roughly equal to the vertical sum of uh.

140 real allocable_, dimension(NIMEM_,NJMEMB_PTR_) :: vhbt !< average y-volume or mass flux determined by the

141 !! barotropic solver [H L2 T-1 ~> m3 s-1 or kg s-1].

142 !! vhbt is roughly equal to vertical sum of vh.

143 real allocable_, dimension(NIMEM_,NJMEM_,NKMEM_) :: pbce !< pbce times eta gives the baroclinic pressure

144 !! anomaly in each layer due to free surface height

145 !! anomalies [L2 H-1 T-2 ~> m s-2 or m4 kg-1 s-2].

146 real allocable_, dimension(NIMEMB_PTR_,NJMEM_) :: du_av_inst !< The barotropic zonal velocity increment

147 !! between filtered and instantaneous velocities

148 !! [L T-1 ~> m s-1]

149 real allocable_, dimension(NIMEM_,NJMEMB_PTR_) :: dv_av_inst !< The barotropic meridional velocity increment

150 !! between filtered and instantaneous velocities

151 !! [L T-1 ~> m s-1]

152

153 real, pointer, dimension(:,:) :: taux_bot => null() !< frictional x-bottom stress from the ocean

154 !! to the seafloor [R L Z T-2 ~> Pa]

155 real, pointer, dimension(:,:) :: tauy_bot => null() !< frictional y-bottom stress from the ocean

156 !! to the seafloor [R L Z T-2 ~> Pa]

157 type(bt_cont_type), pointer :: bt_cont => null() !< A structure with elements that describe the

158 !! effective summed open face areas as a function

159 !! of barotropic flow.

160

161 logical :: bt_adj_corr_mass_src !< If true, recalculates the barotropic mass source after

162 !! predictor step. This should make little difference in the

163 !! deep ocean but appears to help for vanished layers.

164 logical :: split_bottom_stress !< If true, provide the bottom stress

165 !! calculated by the vertical viscosity to the

166 !! barotropic solver.

167 logical :: dtbt_use_bt_cont !< If true, use BT_cont to calculate DTBT.

168 logical :: calculate_sal !< If true, calculate self-attraction and loading.

169 logical :: use_tides !< If true, tidal forcing is enabled.

170 logical :: use_ha !< If true, perform inline harmonic analysis.

171 logical :: remap_aux !< If true, apply ALE remapping to all of the auxiliary 3-D

172 !! variables that are needed to reproduce across restarts,

173 !! similarly to what is done with the primary state variables.

174

175 real :: be !< A nondimensional number from 0.5 to 1 that controls

176 !! the backward weighting of the time stepping scheme [nondim]

177 real :: begw !< A nondimensional number from 0 to 1 that controls

178 !! the extent to which the treatment of gravity waves

179 !! is forward-backward (0) or simulated backward

180 !! Euler (1) [nondim]. 0 is often used.

181 logical :: debug !< If true, write verbose checksums for debugging purposes.

182 logical :: debug_obc !< If true, do additional calls resetting values to help verify the correctness

183 !! of the open boundary condition code.

184 logical :: fpmix = .false. !< If true, applies profiles of momentum flux magnitude and direction.

185 logical :: module_is_initialized = .false. !< Record whether this module has been initialized.

186 logical :: visc_rem_dt_bug = .true. !< If true, recover a bug that uses dt_pred rather than dt for vertvisc_rem

187 !! at the end of predictor.

188

189 !>@{ Diagnostic IDs

190 ! integer :: id_uold = -1, id_vold = -1

191 integer :: id_uh = -1, id_vh = -1

192 integer :: id_umo = -1, id_vmo = -1

193 integer :: id_umo_2d = -1, id_vmo_2d = -1

194 integer :: id_pfu = -1, id_pfv = -1

195 integer :: id_cau = -1, id_cav = -1

196 integer :: id_ueffa = -1, id_veffa = -1

197 ! integer :: id_hf_PFu = -1, id_hf_PFv = -1

198 integer :: id_h_pfu = -1, id_h_pfv = -1

199 integer :: id_hf_pfu_2d = -1, id_hf_pfv_2d = -1

200 integer :: id_intz_pfu_2d = -1, id_intz_pfv_2d = -1

201 integer :: id_pfu_visc_rem = -1, id_pfv_visc_rem = -1

202 ! integer :: id_hf_CAu = -1, id_hf_CAv = -1

203 integer :: id_h_cau = -1, id_h_cav = -1

204 integer :: id_hf_cau_2d = -1, id_hf_cav_2d = -1

205 integer :: id_intz_cau_2d = -1, id_intz_cav_2d = -1

206 integer :: id_cau_visc_rem = -1, id_cav_visc_rem = -1

207 integer :: id_deta_dt = -1

208

209 ! Split scheme only.

210 integer :: id_uav = -1, id_vav = -1

211 integer :: id_u_bt_accel = -1, id_v_bt_accel = -1

212 ! integer :: id_hf_u_BT_accel = -1, id_hf_v_BT_accel = -1

213 integer :: id_h_u_bt_accel = -1, id_h_v_bt_accel = -1

214 integer :: id_hf_u_bt_accel_2d = -1, id_hf_v_bt_accel_2d = -1

215 integer :: id_intz_u_bt_accel_2d = -1, id_intz_v_bt_accel_2d = -1

216 integer :: id_u_bt_accel_visc_rem = -1, id_v_bt_accel_visc_rem = -1

217 !>@}

218

219 type(diag_ctrl), pointer :: diag => null() !< A structure that is used to regulate the

220 !! timing of diagnostic output.

221 type(accel_diag_ptrs), pointer :: adp => null() !< A structure pointing to the various

222 !! accelerations in the momentum equations,

223 !! which can later be used to calculate

224 !! derived diagnostics like energy budgets.

225 type(accel_diag_ptrs), pointer :: ad_pred => null() !< A structure pointing to the various

226 !! predictor step accelerations in the momentum equations,

227 !! which can be used to debug truncations.

228 type(cont_diag_ptrs), pointer :: cdp => null() !< A structure with pointers to various

229 !! terms in the continuity equations,

230 !! which can later be used to calculate

231 !! derived diagnostics like energy budgets.

232

233 ! The remainder of the structure points to child subroutines' control structures.

234 !> A pointer to the horizontal viscosity control structure

235 type(hor_visc_cs) :: hor_visc

236 !> A pointer to the continuity control structure

237 type(continuity_cs) :: continuity_csp

238 !> The CoriolisAdv control structure

239 type(coriolisadv_cs) :: coriolisadv

240 !> A pointer to the PressureForce control structure

241 type(pressureforce_cs) :: pressureforce_csp

242 !> A pointer to a structure containing interface height diffusivities

243 type(vertvisc_cs), pointer :: vertvisc_csp => null()

244 !> A pointer to the set_visc control structure

245 type(set_visc_cs), pointer :: set_visc_csp => null()

246 !> A pointer to the barotropic stepping control structure

247 type(barotropic_cs) :: barotropic_csp

248 !> A pointer to the SAL control structure

249 type(sal_cs) :: sal_csp

250 !> A pointer to the tidal forcing control structure

251 type(tidal_forcing_cs) :: tides_csp

252 !> A pointer to the harmonic analysis control structure

253 type(harmonic_analysis_cs) :: ha_csp

254 !> A pointer to the ALE control structure.

255 type(ale_cs), pointer :: ale_csp => null()

256

257 type(ocean_obc_type), pointer :: obc => null() !< A pointer to an open boundary

258 !! condition type that specifies whether, where, and what open boundary

259 !! conditions are used. If no open BCs are used, this pointer stays

260 !! nullified. Flather OBCs use open boundary_CS as well.

261 !> A pointer to the update_OBC control structure

262 type(update_obc_cs), pointer :: update_obc_csp => null()

263

264 type(group_pass_type) :: pass_eta !< Structure for group halo pass

265 type(group_pass_type) :: pass_visc_rem !< Structure for group halo pass

266 type(group_pass_type) :: pass_uvp !< Structure for group halo pass

267 type(group_pass_type) :: pass_uv_inst !< Structure for group halo pass

268 type(group_pass_type) :: pass_hp_uv !< Structure for group halo pass

269 type(group_pass_type) :: pass_hp_uhvh !< Structure for group halo pass

270 type(group_pass_type) :: pass_h_uv !< Structure for group halo pass

271

272end type mom_dyn_split_rk2b_cs

273

274

275public step_mom_dyn_split_rk2b

276public register_restarts_dyn_split_rk2b

277public initialize_dyn_split_rk2b

278public remap_dyn_split_rk2b_aux_vars

279public end_dyn_split_rk2b

280

281!>@{ CPU time clock IDs

282integer :: id_clock_cor, id_clock_pres, id_clock_vertvisc

283integer :: id_clock_horvisc, id_clock_mom_update

284integer :: id_clock_continuity, id_clock_thick_diff

285integer :: id_clock_btstep, id_clock_btcalc, id_clock_btforce

286integer :: id_clock_pass

287!>@}

288

289contains

290

291!> RK2 splitting for time stepping MOM adiabatic dynamics

292subroutine step_mom_dyn_split_rk2b(u_av, v_av, h, tv, visc, Time_local, dt, forces, &

293 p_surf_begin, p_surf_end, uh, vh, uhtr, vhtr, eta_av, &

294 G, GV, US, CS, calc_dtbt, VarMix, MEKE, thickness_diffuse_CSp, pbv, Waves)

295 type(ocean_grid_type), intent(inout) :: g !< Ocean grid structure

296 type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure

297 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type

298 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &

299 target, intent(inout) :: u_av !< Zonal velocity [L T-1 ~> m s-1]

300 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &

301 target, intent(inout) :: v_av !< Meridional velocity [L T-1 ~> m s-1]

302 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

303 intent(inout) :: h !< Layer thickness [H ~> m or kg m-2]

304 type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic type

305 type(vertvisc_type), intent(inout) :: visc !< Vertical visc, bottom drag, and related

306 type(time_type), intent(in) :: time_local !< Model time at end of time step

307 real, intent(in) :: dt !< Baroclinic dynamics time step [T ~> s]

308 type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces

309 real, dimension(:,:), pointer :: p_surf_begin !< Surface pressure at the start of this dynamic

310 !! time step [R L2 T-2 ~> Pa]

311 real, dimension(:,:), pointer :: p_surf_end !< Surface pressure at the end of this dynamic

312 !! time step [R L2 T-2 ~> Pa]

313 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &

314 target, intent(inout) :: uh !< Zonal volume or mass transport

315 !! [H L2 T-1 ~> m3 s-1 or kg s-1]

316 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &

317 target, intent(inout) :: vh !< Meridional volume or mass transport

318 !! [H L2 T-1 ~> m3 s-1 or kg s-1]

319 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &

320 intent(inout) :: uhtr !< Accumulated zonal volume or mass transport

321 !! since last tracer advection [H L2 ~> m3 or kg]

322 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &

323 intent(inout) :: vhtr !< Accumulated meridional volume or mass transport

324 !! since last tracer advection [H L2 ~> m3 or kg]

325 real, dimension(SZI_(G),SZJ_(G)), intent(out) :: eta_av !< Free surface height or column mass

326 !! averaged over time step [H ~> m or kg m-2]

327 type(mom_dyn_split_rk2b_cs), pointer :: cs !< Module control structure

328 logical, intent(in) :: calc_dtbt !< If true, recalculate the barotropic time step

329 type(varmix_cs), intent(inout) :: varmix !< Variable mixing control structure

330 type(meke_type), intent(inout) :: meke !< MEKE fields

331 type(thickness_diffuse_cs), intent(inout) :: thickness_diffuse_csp !< Pointer to a structure containing

332 !! interface height diffusivities

333 type(porous_barrier_type), intent(in) :: pbv !< porous barrier fractional cell metrics

334 type(wave_parameters_cs), optional, pointer :: waves !< A pointer to a structure containing

335 !! fields related to the surface wave conditions

336

337 ! local variables

338 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: up ! Predicted zonal velocity [L T-1 ~> m s-1].

339 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: vp ! Predicted meridional velocity [L T-1 ~> m s-1].

340 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: hp ! Predicted thickness [H ~> m or kg m-2].

341

342 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: dz ! Distance between the interfaces around a layer [Z ~> m]

343 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: ueffa ! Effective Area of U-Faces [H L ~> m2]

344 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: veffa ! Effective Area of V-Faces [H L ~> m2]

345 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: u_bc_accel ! The summed zonal baroclinic accelerations

346 ! of each layer calculated by the non-barotropic

347 ! part of the model [L T-2 ~> m s-2]

348 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: v_bc_accel ! The summed meridional baroclinic accelerations

349 ! of each layer calculated by the non-barotropic

350 ! part of the model [L T-2 ~> m s-2]

351

352 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), target :: u_inst ! Instantaneous zonal velocity [L T-1 ~> m s-1].

353 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), target :: v_inst ! Instantaneous meridional velocity [L T-1 ~> m s-1].

354 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_av ! The average of the layer thicknesses at the beginning

355 ! and end of a time step [H ~> m or kg m-2]

356

357 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), target :: uh_in ! The zonal mass transports that would be

358 ! obtained using the initial velocities [H L2 T-1 ~> m3 s-1 or kg s-1]

359 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), target :: vh_in ! The meridional mass transports that would be

360 ! obtained using the initial velocities [H L2 T-1 ~> m3 s-1 or kg s-1]

361

362 real, dimension(SZI_(G),SZJ_(G)) :: eta_pred ! The predictor value of the free surface height

363 ! or column mass [H ~> m or kg m-2]

364 real, dimension(SZI_(G),SZJ_(G)) :: spv_avg ! The column averaged specific volume [R-1 ~> m3 kg-1]

365 real, dimension(SZI_(G),SZJ_(G)) :: deta_dt ! A diagnostic of the time derivative of the free surface

366 ! height or column mass [H T-1 ~> m s-1 or kg m-2 s-1]

367

368 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: u_old_rad_obc ! The starting zonal velocities, which are

369 ! saved for use in the Flather open boundary condition code [L T-1 ~> m s-1]

370 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: v_old_rad_obc ! The starting meridional velocities, which are

371 ! saved for use in the Flather open boundary condition code [L T-1 ~> m s-1]

372

373 ! GMM, TODO: make these allocatable?

374 ! real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: uold ! u-velocity before vert_visc is applied, for fpmix

375 ! ! [L T-1 ~> m s-1]

376 ! real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: vold ! v-velocity before vert_visc is applied, for fpmix

377 ! ! [L T-1 ~> m s-1]

378 real :: pres_to_eta ! A factor that converts pressures to the units of eta

379 ! [H T2 R-1 L-2 ~> m Pa-1 or kg m-2 Pa-1]

380 real, pointer, dimension(:,:) :: &

381 p_surf => null(), & ! A pointer to the surface pressure [R L2 T-2 ~> Pa]

382 eta_pf_start => null(), & ! The value of eta that corresponds to the starting pressure

383 ! for the barotropic solver [H ~> m or kg m-2]

384 taux_bot => null(), & ! A pointer to the zonal bottom stress in some cases [R L Z T-2 ~> Pa]

385 tauy_bot => null(), & ! A pointer to the meridional bottom stress in some cases [R L Z T-2 ~> Pa]

386 ! This pointer is just used as shorthand for CS%eta.

387 eta => null() ! A pointer to the instantaneous free surface height (in Boussinesq

388 ! mode) or column mass anomaly (in non-Boussinesq mode) [H ~> m or kg m-2]

389

390 real, pointer, dimension(:,:,:) :: &

391 ! These pointers are used to alter which fields are passed to btstep with various options:

392 u_ptr => null(), & ! A pointer to a zonal velocity [L T-1 ~> m s-1]

393 v_ptr => null(), & ! A pointer to a meridional velocity [L T-1 ~> m s-1]

394 uh_ptr => null(), & ! A pointer to a zonal volume or mass transport [H L2 T-1 ~> m3 s-1 or kg s-1]

395 vh_ptr => null() ! A pointer to a meridional volume or mass transport [H L2 T-1 ~> m3 s-1 or kg s-1]

396

397 ! real, dimension(SZI_(G),SZJ_(G)) :: hbl ! Boundary layer depth from Cvmix [H ~> m or kg m-2]

398 real :: dt_pred ! The time step for the predictor part of the baroclinic time stepping [T ~> s].

399 real :: idt_bc ! Inverse of the baroclinic timestep [T-1 ~> s-1]

400 logical :: dyn_p_surf

401 logical :: debug_redundant ! If true, check redundant values on PE boundaries when debugging

402 logical :: bt_cont_bt_thick ! If true, use the BT_cont_type to estimate the

403 ! relative weightings of the layers in calculating

404 ! the barotropic accelerations.

405 !---For group halo pass

406 logical :: showcalltree, sym

407

408 integer :: i, j, k, is, ie, js, je, isq, ieq, jsq, jeq, nz

409 integer :: cont_stencil, obc_stencil

410 integer :: cor_stencil

411

412 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke

413 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB

414

415 eta => cs%eta

416

417 idt_bc = 1.0 / dt

418

419 sym = g%Domain%symmetric ! switch to include symmetric domain in checksums

420

421 showcalltree = calltree_showquery()

422 if (showcalltree) call calltree_enter("step_MOM_dyn_split_RK2b(), MOM_dynamics_split_RK2b.F90")

423

424 ! Fill in some halo points for arrays that will have halo updates.

425 hp(:,:,:) = h(:,:,:)

426 up(:,:,:) = 0.0 ; vp(:,:,:) = 0.0 ; u_inst(:,:,:) = 0.0 ; v_inst(:,:,:) = 0.0

427

428 ! Update CFL truncation value as function of time

429 call updatecfltruncationvalue(time_local, cs%vertvisc_CSp, us)

430

431 if (cs%debug) then

432 call query_debugging_checks(do_redundant=debug_redundant)

433 call mom_state_chksum("Start predictor ", u_av, v_av, h, uh, vh, g, gv, us, symmetric=sym)

434 if (debug_redundant) then

435 call check_redundant("Start predictor u ", u_av, v_av, g, unscale=us%L_T_to_m_s)

436 call check_redundant("Start predictor uh ", uh, vh, g, unscale=gv%H_to_MKS*us%L_to_m**2*us%s_to_T)

437 endif

438 endif

439

440 dyn_p_surf = associated(p_surf_begin) .and. associated(p_surf_end)

441 if (dyn_p_surf) then

442 p_surf => p_surf_end

443 call safe_alloc_ptr(eta_pf_start,g%isd,g%ied,g%jsd,g%jed)

444 eta_pf_start(:,:) = 0.0

445 else

446 p_surf => forces%p_surf

447 endif

448

449 if (associated(cs%OBC)) then

450 if (cs%debug_OBC) call open_boundary_test_extern_h(g, gv, cs%OBC, h)

451

452 ! Update OBC ramp value as function of time

453 call update_obc_ramp(time_local, cs%OBC, us)

454

455 do k=1,nz ; do j=g%jsd,g%jed ; do i=g%IsdB,g%IedB

456 u_old_rad_obc(i,j,k) = u_av(i,j,k)

457 enddo ; enddo ; enddo

458 do k=1,nz ; do j=g%JsdB,g%JedB ; do i=g%isd,g%ied

459 v_old_rad_obc(i,j,k) = v_av(i,j,k)

460 enddo ; enddo ; enddo

461 endif

462

463 bt_cont_bt_thick = .false.

464 if (associated(cs%BT_cont)) bt_cont_bt_thick = &

465 (allocated(cs%BT_cont%h_u) .and. allocated(cs%BT_cont%h_v))

466

467 if (cs%split_bottom_stress) then

468 taux_bot => cs%taux_bot ; tauy_bot => cs%tauy_bot

469 endif

470

471 !--- begin set up for group halo pass

472

473 cont_stencil = continuity_stencil(cs%continuity_CSp)

474 obc_stencil = 2

475 cor_stencil = coriolisadv_stencil(cs%CoriolisAdv)

476 if (associated(cs%OBC)) then

477 if (cs%OBC%oblique_BCs_exist_globally) obc_stencil = 3

478 endif

479 call cpu_clock_begin(id_clock_pass)

480 call create_group_pass(cs%pass_eta, eta, g%Domain, halo=1)

481 call create_group_pass(cs%pass_visc_rem, cs%visc_rem_u, cs%visc_rem_v, g%Domain, &

482 to_all+scalar_pair, cgrid_ne, halo=max(1,cont_stencil))

483 call create_group_pass(cs%pass_uvp, up, vp, g%Domain, halo=max(1,cont_stencil))

484 call create_group_pass(cs%pass_uv_inst, u_inst, v_inst, g%Domain, halo=max(2,cont_stencil))

485

486 call create_group_pass(cs%pass_hp_uv, hp, g%Domain, halo=cor_stencil)

487 call create_group_pass(cs%pass_hp_uv, u_av, v_av, g%Domain, halo=max(cor_stencil,obc_stencil))

488 call create_group_pass(cs%pass_hp_uv, uh(:,:,:), vh(:,:,:), g%Domain, halo=max(cor_stencil,obc_stencil))

489

490 call create_group_pass(cs%pass_hp_uhvh, hp, g%Domain, halo=cor_stencil)

491 call create_group_pass(cs%pass_hp_uhvh, uh(:,:,:), vh(:,:,:), g%Domain, halo=max(cor_stencil,obc_stencil))

492 if (cor_stencil > 2) then

493 call create_group_pass(cs%pass_hp_uhvh, u_av, v_av, g%Domain, halo=max(cor_stencil,obc_stencil))

494 call create_group_pass(cs%pass_hp_uhvh, h, g%Domain, halo=cor_stencil)

495 endif

496

497 call create_group_pass(cs%pass_h_uv, h, g%Domain, halo=max(cor_stencil,cont_stencil))

498 call create_group_pass(cs%pass_h_uv, u_av, v_av, g%Domain, halo=max(cor_stencil,obc_stencil))

499 call create_group_pass(cs%pass_h_uv, uh(:,:,:), vh(:,:,:), g%Domain, halo=max(cor_stencil,obc_stencil))

500

501 call cpu_clock_end(id_clock_pass)

502 !--- end set up for group halo pass

503

504 ! This calculates the transports and averaged thicknesses that will be used for the

505 ! predictor version of the Coriolis scheme.

506 call cpu_clock_begin(id_clock_continuity)

507 call continuity(u_av, v_av, h, hp, uh, vh, dt, g, gv, us, cs%continuity_CSp, cs%OBC, pbv)

508 call cpu_clock_end(id_clock_continuity)

509

510 if (g%nonblocking_updates) &

511 call start_group_pass(cs%pass_hp_uhvh, g%Domain, clock=id_clock_pass)

512

513! PFu = d/dx M(h,T,S)

514! pbce = dM/deta

515 if (cs%begw == 0.0) call enable_averages(dt, time_local, cs%diag)

516 call cpu_clock_begin(id_clock_pres)

517 call pressureforce(h, tv, cs%PFu, cs%PFv, g, gv, us, cs%PressureForce_CSp, &

518 cs%ALE_CSp, cs%ADp, p_surf, cs%pbce, cs%eta_PF)

519 if (dyn_p_surf) then

520 pres_to_eta = 1.0 / (gv%g_Earth * gv%H_to_RZ)

521 !$OMP parallel do default(shared)

522 do j=jsq,jeq+1 ; do i=isq,ieq+1

523 eta_pf_start(i,j) = cs%eta_PF(i,j) - pres_to_eta * (p_surf_begin(i,j) - p_surf_end(i,j))

524 enddo ; enddo

525 endif

526 ! Stokes shear force contribution to pressure gradient

527 if (present(waves)) then ; if (associated(waves)) then ; if (waves%Stokes_PGF) then

528 call thickness_to_dz(h, tv, dz, g, gv, us, halo_size=1)

529 call stokes_pgf(g, gv, us, dz, u_av, v_av, cs%PFu_Stokes, cs%PFv_Stokes, waves)

530

531 ! We are adding Stokes_PGF to hydrostatic PGF here. The diag PFu/PFv

532 ! will therefore report the sum total PGF and we avoid other

533 ! modifications in the code. The PFu_Stokes is output within the waves routines.

534 if (.not.waves%Passive_Stokes_PGF) then

535 do k=1,nz ; do j=js,je ; do i=isq,ieq

536 cs%PFu(i,j,k) = cs%PFu(i,j,k) + cs%PFu_Stokes(i,j,k)

537 enddo ; enddo ; enddo

538 do k=1,nz ; do j=jsq,jeq ; do i=is,ie

539 cs%PFv(i,j,k) = cs%PFv(i,j,k) + cs%PFv_Stokes(i,j,k)

540 enddo ; enddo ; enddo

541 endif

542 endif ; endif ; endif

543 call cpu_clock_end(id_clock_pres)

544 call disable_averaging(cs%diag)

545 if (showcalltree) call calltree_waypoint("done with PressureForce (step_MOM_dyn_split_RK2b)")

546

547 if (associated(cs%OBC)) then ; if (cs%OBC%update_OBC) then

548 call update_obc_data(cs%OBC, g, gv, us, tv, h, cs%update_OBC_CSp, time_local)

549 endif ; endif

550 if (associated(cs%OBC) .and. cs%debug_OBC) &

551 call open_boundary_zero_normal_flow(cs%OBC, g, gv, cs%PFu, cs%PFv)

552

553 if (g%nonblocking_updates) then

554 call complete_group_pass(cs%pass_hp_uhvh, g%Domain, clock=id_clock_pass)

555 call start_group_pass(cs%pass_eta, g%Domain, clock=id_clock_pass)

556 else

557 call do_group_pass(cs%pass_hp_uhvh, g%Domain, clock=id_clock_pass)

558 endif

559

560 !$OMP parallel do default(shared)

561 do k=1,nz ; do j=js-cor_stencil,je+cor_stencil ; do i=is-cor_stencil,ie+cor_stencil

562 h_av(i,j,k) = 0.5*(h(i,j,k) + hp(i,j,k))

563 enddo ; enddo ; enddo

564

565 ! Calculate a predictor-step estimate of the Coriolis and momentum advection terms

566 ! and horizontal viscous accelerations.

567

568! CAu = -(f+zeta_av)/h_av vh + d/dx KE_av

569 call cpu_clock_begin(id_clock_cor)

570 call coradcalc(u_av, v_av, h_av, uh, vh, cs%CAu_pred, cs%CAv_pred, cs%OBC, cs%AD_pred, &

571 g, gv, us, cs%CoriolisAdv, pbv, waves=waves)

572 call cpu_clock_end(id_clock_cor)

573 if (showcalltree) call calltree_waypoint("done with predictor CorAdCalc (step_MOM_dyn_split_RK2b)")

574

575! diffu = horizontal viscosity terms (u_av)

576 call cpu_clock_begin(id_clock_horvisc)

577 call horizontal_viscosity(u_av, v_av, h_av, uh, vh, cs%diffu, cs%diffv, meke, varmix, g, gv, us, cs%hor_visc, &

578 tv, dt, obc=cs%OBC, bt=cs%barotropic_CSp, td=thickness_diffuse_csp, adp=cs%AD_pred)

579 call cpu_clock_end(id_clock_horvisc)

580 if (showcalltree) call calltree_waypoint("done with predictor horizontal_viscosity (step_MOM_dyn_split_RK2b)")

581

582! u_bc_accel = CAu + PFu + diffu(u[n-1])

583 call cpu_clock_begin(id_clock_btforce)

584 !$OMP parallel do default(shared)

585 do k=1,nz

586 do j=js,je ; do i=isq,ieq

587 u_bc_accel(i,j,k) = (cs%CAu_pred(i,j,k) + cs%PFu(i,j,k)) + cs%diffu(i,j,k)

588 enddo ; enddo

589 do j=jsq,jeq ; do i=is,ie

590 v_bc_accel(i,j,k) = (cs%CAv_pred(i,j,k) + cs%PFv(i,j,k)) + cs%diffv(i,j,k)

591 enddo ; enddo

592 enddo

593 if (associated(cs%OBC)) then

594 call open_boundary_zero_normal_flow(cs%OBC, g, gv, u_bc_accel, v_bc_accel)

595 endif

596 call cpu_clock_end(id_clock_btforce)

597

598 if (cs%debug) then

599 call mom_accel_chksum("pre-btstep accel", cs%CAu_pred, cs%CAv_pred, cs%PFu, cs%PFv, &

600 cs%diffu, cs%diffv, g, gv, us, cs%pbce, u_bc_accel, v_bc_accel, &

601 symmetric=sym)

602 if (debug_redundant) then

603 call check_redundant("pre-btstep CS%CA ", cs%CAu_pred, cs%CAv_pred, g, unscale=us%L_T2_to_m_s2)

604 call check_redundant("pre-btstep CS%PF ", cs%PFu, cs%PFv, g, unscale=us%L_T2_to_m_s2)

605 call check_redundant("pre-btstep CS%diff ", cs%diffu, cs%diffv, g, unscale=us%L_T2_to_m_s2)

606 call check_redundant("pre-btstep u_bc_accel ", u_bc_accel, v_bc_accel, g, unscale=us%L_T2_to_m_s2)

607 endif

608 endif

609

610 call cpu_clock_begin(id_clock_vertvisc)

611 !$OMP parallel do default(shared)

612 do k=1,nz

613 do j=js,je ; do i=isq,ieq

614 up(i,j,k) = g%mask2dCu(i,j) * (u_av(i,j,k) + dt * u_bc_accel(i,j,k))

615 enddo ; enddo

616 do j=jsq,jeq ; do i=is,ie

617 vp(i,j,k) = g%mask2dCv(i,j) * (v_av(i,j,k) + dt * v_bc_accel(i,j,k))

618 enddo ; enddo

619 enddo

620

621 call enable_averages(dt, time_local, cs%diag)

622 call set_viscous_ml(u_av, v_av, h, tv, forces, visc, dt, g, gv, us, cs%set_visc_CSp)

623 call disable_averaging(cs%diag)

624

625 if (cs%debug) then

626 call uvchksum("before vertvisc: up", up, vp, g%HI, haloshift=0, symmetric=sym, unscale=us%L_T_to_m_s)

627 endif

628 call thickness_to_dz(h, tv, dz, g, gv, us, halo_size=1)

629 call vertvisc_coef(up, vp, h, dz, forces, visc, tv, dt, g, gv, us, cs%vertvisc_CSp, cs%OBC, varmix)

630 call vertvisc_remnant(visc, cs%visc_rem_u, cs%visc_rem_v, dt, g, gv, us, cs%vertvisc_CSp)

631 call cpu_clock_end(id_clock_vertvisc)

632 if (showcalltree) call calltree_waypoint("done with vertvisc_coef (step_MOM_dyn_split_RK2b)")

633

634

635 call cpu_clock_begin(id_clock_pass)

636 if (g%nonblocking_updates) then

637 call complete_group_pass(cs%pass_eta, g%Domain)

638 call start_group_pass(cs%pass_visc_rem, g%Domain)

639 else

640 call do_group_pass(cs%pass_eta, g%Domain)

641 call do_group_pass(cs%pass_visc_rem, g%Domain)

642 endif

643 call cpu_clock_end(id_clock_pass)

644

645 call cpu_clock_begin(id_clock_btcalc)

646 ! Calculate the relative layer weights for determining barotropic quantities.

647 if (.not.bt_cont_bt_thick) &

648 call btcalc(h, g, gv, cs%barotropic_CSp, obc=cs%OBC)

649 call bt_mass_source(h, eta, .true., g, gv, cs%barotropic_CSp)

650

651 spv_avg(:,:) = 0.0

652 if ((.not.gv%Boussinesq) .and. associated(cs%OBC)) then

653 ! Determine the column average specific volume if it is needed due to the

654 ! use of Flather open boundary conditions in non-Boussinesq mode.

655 if (open_boundary_query(cs%OBC, apply_flather_obc=.true.)) &

656 call find_col_avg_spv(h, spv_avg, tv, g, gv, us)

657 endif

658 call cpu_clock_end(id_clock_btcalc)

659

660 if (g%nonblocking_updates) &

661 call complete_group_pass(cs%pass_visc_rem, g%Domain, clock=id_clock_pass)

662

663 ! Reconstruct u_inst and v_inst from u_av and v_av.

664 call cpu_clock_begin(id_clock_mom_update)

665 do k=1,nz ; do j=js,je ; do i=isq,ieq

666 u_inst(i,j,k) = u_av(i,j,k) - cs%du_av_inst(i,j) * cs%visc_rem_u(i,j,k)

667 enddo ; enddo ; enddo

668 do k=1,nz ; do j=jsq,jeq ; do i=is,ie

669 v_inst(i,j,k) = v_av(i,j,k) - cs%dv_av_inst(i,j) * cs%visc_rem_v(i,j,k)

670 enddo ; enddo ; enddo

671 call cpu_clock_end(id_clock_mom_update)

672 call do_group_pass(cs%pass_uv_inst, g%Domain, clock=id_clock_pass)

673

674 if (associated(cs%OBC)) &

675 call copy_thickness_reservoirs(cs%OBC, g, gv)

676

677! u_accel_bt = layer accelerations due to barotropic solver

678 call cpu_clock_begin(id_clock_continuity)

679 call continuity(u_inst, v_inst, h, hp, uh_in, vh_in, dt, g, gv, us, cs%continuity_CSp, cs%OBC, pbv, &

680 visc_rem_u=cs%visc_rem_u, visc_rem_v=cs%visc_rem_v, bt_cont=cs%BT_cont)

681 call cpu_clock_end(id_clock_continuity)

682 if (bt_cont_bt_thick) then

683 call btcalc(h, g, gv, cs%barotropic_CSp, cs%BT_cont%h_u, cs%BT_cont%h_v, &

684 obc=cs%OBC)

685 endif

686 if (showcalltree) call calltree_waypoint("done with continuity[BT_cont] (step_MOM_dyn_split_RK2b)")

687

688 uh_ptr => uh_in ; vh_ptr => vh_in ; u_ptr => u_inst ; v_ptr => v_inst

689

690 call cpu_clock_begin(id_clock_btstep)

691 if (calc_dtbt) then

692 if (cs%dtbt_use_bt_cont .and. associated(cs%BT_cont)) then

693 call set_dtbt(g, gv, us, cs%barotropic_CSp, cs%pbce, bt_cont=cs%BT_cont)

694 else

695 ! In the following call, eta is only used when NONLINEAR_BT_CONTINUITY is True. Otherwise, dtbt is effectively

696 ! calculated with eta=0. Note that NONLINEAR_BT_CONTINUITY is False if BT_CONT is used, which is the default.

697 call set_dtbt(g, gv, us, cs%barotropic_CSp, cs%pbce, eta=eta)

698 endif

699 endif

700 if (showcalltree) call calltree_enter("btstep(), MOM_barotropic.F90")

701 ! This is the predictor step call to btstep.

702 ! The CS%ADp argument here stores the weights for certain integrated diagnostics.

703 call btstep(u_inst, v_inst, eta, dt, u_bc_accel, v_bc_accel, forces, cs%pbce, cs%eta_PF, u_av, v_av, &

704 cs%u_accel_bt, cs%v_accel_bt, eta_pred, cs%uhbt, cs%vhbt, g, gv, us, &

705 cs%barotropic_CSp, cs%visc_rem_u, cs%visc_rem_v, spv_avg, cs%ADp, cs%OBC, cs%BT_cont, &

706 eta_pf_start, taux_bot, tauy_bot, uh_ptr, vh_ptr, u_ptr, v_ptr)

707 if (showcalltree) call calltree_leave("btstep()")

708 call cpu_clock_end(id_clock_btstep)

709

710 dt_pred = dt * cs%be

711 call cpu_clock_begin(id_clock_mom_update)

712

713! up = u + dt_pred*( u_bc_accel + u_accel_bt )

714 !$OMP parallel do default(shared)

715 do k=1,nz

716 do j=jsq,jeq ; do i=is,ie

717 vp(i,j,k) = g%mask2dCv(i,j) * (v_inst(i,j,k) + dt_pred * &

718 (v_bc_accel(i,j,k) + cs%v_accel_bt(i,j,k)))

719 enddo ; enddo

720 do j=js,je ; do i=isq,ieq

721 up(i,j,k) = g%mask2dCu(i,j) * (u_inst(i,j,k) + dt_pred * &

722 (u_bc_accel(i,j,k) + cs%u_accel_bt(i,j,k)))

723 enddo ; enddo

724 enddo

725 call cpu_clock_end(id_clock_mom_update)

726

727 if (cs%debug) then

728 call uvchksum("Predictor 1 [uv]", up, vp, g%HI, haloshift=0, symmetric=sym, unscale=us%L_T_to_m_s)

729 call hchksum(h, "Predictor 1 h", g%HI, haloshift=1, unscale=gv%H_to_MKS)

730 call uvchksum("Predictor 1 [uv]h", uh, vh, g%HI,haloshift=2, &

731 symmetric=sym, unscale=gv%H_to_MKS*us%L_to_m**2*us%s_to_T)

732! call MOM_state_chksum("Predictor 1", up, vp, h, uh, vh, G, GV, US, haloshift=1)

733 call mom_accel_chksum("Predictor accel", cs%CAu_pred, cs%CAv_pred, cs%PFu, cs%PFv, &

734 cs%diffu, cs%diffv, g, gv, us, cs%pbce, cs%u_accel_bt, cs%v_accel_bt, symmetric=sym)

735 call mom_state_chksum("Predictor 1 init", u_inst, v_inst, h, uh, vh, g, gv, us, haloshift=1, &

736 symmetric=sym)

737 if (debug_redundant) then

738 call check_redundant("Predictor 1 up", up, vp, g, unscale=us%L_T_to_m_s)

739 call check_redundant("Predictor 1 uh", uh, vh, g, unscale=gv%H_to_MKS*us%L_to_m**2*us%s_to_T)

740 endif

741 endif

742

743! up <- up + dt_pred d/dz visc d/dz up

744! u_av <- u_av + dt_pred d/dz visc d/dz u_av

745 call cpu_clock_begin(id_clock_vertvisc)

746 if (cs%debug) then

747 call uvchksum("0 before vertvisc: [uv]p", up, vp, g%HI,haloshift=0, symmetric=sym, unscale=us%L_T_to_m_s)

748 endif

749

750 ! if (CS%fpmix) then

751 ! uold(:,:,:) = 0.0

752 ! vold(:,:,:) = 0.0

753 ! do k=1,nz ; do j=js,je ; do I=Isq,Ieq

754 ! uold(I,j,k) = up(I,j,k)

755 ! enddo ; enddo ; enddo

756 ! do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie

757 ! vold(i,J,k) = vp(i,J,k)

758 ! enddo ; enddo ; enddo

759 ! endif

760

761 call thickness_to_dz(h, tv, dz, g, gv, us, halo_size=1)

762 call vertvisc_coef(up, vp, h, dz, forces, visc, tv, dt_pred, g, gv, us, cs%vertvisc_CSp, &

763 cs%OBC, varmix)

764 call vertvisc(up, vp, h, forces, visc, dt_pred, cs%OBC, cs%AD_pred, cs%CDp, g, &

765 gv, us, cs%vertvisc_CSp, cs%taux_bot, cs%tauy_bot, waves=waves)

766

767 ! if (CS%fpmix) then

768 ! hbl(:,:) = 0.0

769 ! if (associated(visc%h_ML)) hbl(:,:) = visc%h_ML(:,:)

770 ! call vertFPmix(up, vp, uold, vold, hbl, h, forces, &

771 ! dt_pred, G, GV, US, CS%vertvisc_CSp, CS%OBC)

772 ! call vertvisc(up, vp, h, forces, visc, dt_pred, CS%OBC, CS%ADp, CS%CDp, G, &

773 ! GV, US, CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot, waves=waves)

774 ! endif

775

776 if (showcalltree) call calltree_waypoint("done with vertvisc (step_MOM_dyn_split_RK2b)")

777 if (g%nonblocking_updates) then

778 call cpu_clock_end(id_clock_vertvisc)

779 call start_group_pass(cs%pass_uvp, g%Domain, clock=id_clock_pass)

780 call cpu_clock_begin(id_clock_vertvisc)

781 endif

782 if (cs%visc_rem_dt_bug) then

783 call vertvisc_remnant(visc, cs%visc_rem_u, cs%visc_rem_v, dt_pred, g, gv, us, cs%vertvisc_CSp)

784 else

785 call vertvisc_remnant(visc, cs%visc_rem_u, cs%visc_rem_v, dt, g, gv, us, cs%vertvisc_CSp)

786 endif

787 call cpu_clock_end(id_clock_vertvisc)

788

789 call do_group_pass(cs%pass_visc_rem, g%Domain, clock=id_clock_pass)

790 if (g%nonblocking_updates) then

791 call complete_group_pass(cs%pass_uvp, g%Domain, clock=id_clock_pass)

792 else

793 call do_group_pass(cs%pass_uvp, g%Domain, clock=id_clock_pass)

794 endif

795

796 ! uh = u_av * h

797 ! hp = h + dt * div . uh

798 call cpu_clock_begin(id_clock_continuity)

799 call continuity(up, vp, h, hp, uh, vh, dt, g, gv, us, cs%continuity_CSp, cs%OBC, pbv, &

800 cs%uhbt, cs%vhbt, cs%visc_rem_u, cs%visc_rem_v, &

801 u_av, v_av, bt_cont=cs%BT_cont)

802 call cpu_clock_end(id_clock_continuity)

803 if (showcalltree) call calltree_waypoint("done with continuity (step_MOM_dyn_split_RK2b)")

804

805 call do_group_pass(cs%pass_hp_uv, g%Domain, clock=id_clock_pass)

806

807 if (associated(cs%OBC)) then

808

809 if (cs%debug) &

810 call uvchksum("Pre OBC avg [uv]", u_av, v_av, g%HI, haloshift=1, symmetric=sym, unscale=us%L_T_to_m_s)

811

812 call radiation_open_bdry_conds(cs%OBC, u_av, u_old_rad_obc, v_av, v_old_rad_obc, g, gv, us, dt_pred)

813

814 if (cs%debug) &

815 call uvchksum("Post OBC avg [uv]", u_av, v_av, g%HI, haloshift=1, symmetric=sym, unscale=us%L_T_to_m_s)

816 endif

817

818 ! h_av = (h + hp)/2

819 !$OMP parallel do default(shared)

820 do k=1,nz ; do j=js-cor_stencil,je+cor_stencil ; do i=is-cor_stencil,ie+cor_stencil

821 h_av(i,j,k) = 0.5*(h(i,j,k) + hp(i,j,k))

822 enddo ; enddo ; enddo

823

824 ! The correction phase of the time step starts here.

825 call enable_averages(dt, time_local, cs%diag)

826

827 ! Calculate a revised estimate of the free-surface height correction to be

828 ! used in the next call to btstep. This call is at this point so that

829 ! hp can be changed if CS%begw /= 0.

830 ! eta_cor = ... (hidden inside CS%barotropic_CSp)

831 if (cs%BT_adj_corr_mass_src) then

832 call cpu_clock_begin(id_clock_btcalc)

833 call bt_mass_source(hp, eta_pred, .false., g, gv, cs%barotropic_CSp)

834 call cpu_clock_end(id_clock_btcalc)

835 endif

836

837 if (cs%begw /= 0.0) then

838 ! hp <- (1-begw)*h_in + begw*hp

839 ! Back up hp to the value it would have had after a time-step of

840 ! begw*dt. hp is not used again until recalculated by continuity.

841 !$OMP parallel do default(shared)

842 do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1

843 hp(i,j,k) = (1.0-cs%begw)*h(i,j,k) + cs%begw*hp(i,j,k)

844 enddo ; enddo ; enddo

845

846 ! PFu = d/dx M(hp,T,S)

847 ! pbce = dM/deta

848 call cpu_clock_begin(id_clock_pres)

849 call pressureforce(hp, tv, cs%PFu, cs%PFv, g, gv, us, cs%PressureForce_CSp, &

850 cs%ALE_CSp, cs%ADp, p_surf, cs%pbce, cs%eta_PF)

851 ! Stokes shear force contribution to pressure gradient

852 if (present(waves)) then ; if (associated(waves)) then ; if (waves%Stokes_PGF) then

853 call thickness_to_dz(h, tv, dz, g, gv, us, halo_size=1)

854 call stokes_pgf(g, gv, us, dz, u_av, v_av, cs%PFu_Stokes, cs%PFv_Stokes, waves)

855 if (.not.waves%Passive_Stokes_PGF) then

856 do k=1,nz ; do j=js,je ; do i=isq,ieq

857 cs%PFu(i,j,k) = cs%PFu(i,j,k) + cs%PFu_Stokes(i,j,k)

858 enddo ; enddo ; enddo

859 do k=1,nz ; do j=jsq,jeq ; do i=is,ie

860 cs%PFv(i,j,k) = cs%PFv(i,j,k) + cs%PFv_Stokes(i,j,k)

861 enddo ; enddo ; enddo

862 endif

863 endif ; endif ; endif

864 call cpu_clock_end(id_clock_pres)

865 if (showcalltree) call calltree_waypoint("done with PressureForce[hp=(1-b).h+b.h] (step_MOM_dyn_split_RK2b)")

866 endif

867

868 if (bt_cont_bt_thick) then

869 call btcalc(h, g, gv, cs%barotropic_CSp, cs%BT_cont%h_u, cs%BT_cont%h_v, &

870 obc=cs%OBC)

871 if (showcalltree) call calltree_waypoint("done with btcalc[BT_cont_BT_thick] (step_MOM_dyn_split_RK2b)")

872 endif

873

874 if (cs%debug) then

875 call mom_state_chksum("Predictor ", up, vp, hp, uh, vh, g, gv, us, symmetric=sym)

876 call uvchksum("Predictor avg [uv]", u_av, v_av, g%HI, haloshift=1, symmetric=sym, unscale=us%L_T_to_m_s)

877 call hchksum(h_av, "Predictor avg h", g%HI, haloshift=2, unscale=gv%H_to_MKS)

878 ! call MOM_state_chksum("Predictor avg ", u_av, v_av, h_av, uh, vh, G, GV, US)

879 if (debug_redundant) then

880 call check_redundant("Predictor up ", up, vp, g, unscale=us%L_T_to_m_s)

881 call check_redundant("Predictor uh ", uh, vh, g, unscale=gv%H_to_MKS*us%L_to_m**2*us%s_to_T)

882 endif

883 endif

884

885! diffu = horizontal viscosity terms (u_av)

886 call cpu_clock_begin(id_clock_horvisc)

887 call horizontal_viscosity(u_av, v_av, h_av, uh, vh, cs%diffu, cs%diffv, meke, varmix, g, gv, us, cs%hor_visc, &

888 tv, dt, obc=cs%OBC, bt=cs%barotropic_CSp, td=thickness_diffuse_csp, adp=cs%ADp)

889 call cpu_clock_end(id_clock_horvisc)

890 if (showcalltree) call calltree_waypoint("done with horizontal_viscosity (step_MOM_dyn_split_RK2b)")

891

892! CAu = -(f+zeta_av)/h_av vh + d/dx KE_av

893 call cpu_clock_begin(id_clock_cor)

894 call coradcalc(u_av, v_av, h_av, uh, vh, cs%CAu, cs%CAv, cs%OBC, cs%ADp, &

895 g, gv, us, cs%CoriolisAdv, pbv, waves=waves)

896 call cpu_clock_end(id_clock_cor)

897 if (showcalltree) call calltree_waypoint("done with CorAdCalc (step_MOM_dyn_split_RK2b)")

898

899! Calculate the momentum forcing terms for the barotropic equations.

900

901! u_bc_accel = CAu + PFu + diffu(u[n-1])

902 call cpu_clock_begin(id_clock_btforce)

903 !$OMP parallel do default(shared)

904 do k=1,nz

905 do j=js,je ; do i=isq,ieq

906 u_bc_accel(i,j,k) = (cs%Cau(i,j,k) + cs%PFu(i,j,k)) + cs%diffu(i,j,k)

907 enddo ; enddo

908 do j=jsq,jeq ; do i=is,ie

909 v_bc_accel(i,j,k) = (cs%Cav(i,j,k) + cs%PFv(i,j,k)) + cs%diffv(i,j,k)

910 enddo ; enddo

911 enddo

912 if (associated(cs%OBC)) then

913 call open_boundary_zero_normal_flow(cs%OBC, g, gv, u_bc_accel, v_bc_accel)

914 endif

915 call cpu_clock_end(id_clock_btforce)

916

917 if (cs%debug) then

918 call mom_accel_chksum("corr pre-btstep accel", cs%CAu, cs%CAv, cs%PFu, cs%PFv, &

919 cs%diffu, cs%diffv, g, gv, us, cs%pbce, u_bc_accel, v_bc_accel, &

920 symmetric=sym)

921 if (debug_redundant) then

922 call check_redundant("corr pre-btstep CS%CA ", cs%CAu, cs%CAv, g, unscale=us%L_T2_to_m_s2)

923 call check_redundant("corr pre-btstep CS%PF ", cs%PFu, cs%PFv, g, unscale=us%L_T2_to_m_s2)

924 call check_redundant("corr pre-btstep CS%diff ", cs%diffu, cs%diffv, g, unscale=us%L_T2_to_m_s2)

925 call check_redundant("corr pre-btstep u_bc_accel ", u_bc_accel, v_bc_accel, g, unscale=us%L_T2_to_m_s2)

926 endif

927 endif

928

929 ! u_accel_bt = layer accelerations due to barotropic solver

930 ! pbce = dM/deta

931 call cpu_clock_begin(id_clock_btstep)

932

933 uh_ptr => uh ; vh_ptr => vh ; u_ptr => u_av ; v_ptr => v_av

934

935 if (showcalltree) call calltree_enter("btstep(), MOM_barotropic.F90")

936 ! This is the corrector step call to btstep.

937 call btstep(u_inst, v_inst, eta, dt, u_bc_accel, v_bc_accel, forces, cs%pbce, cs%eta_PF, u_av, v_av, &

938 cs%u_accel_bt, cs%v_accel_bt, eta_pred, cs%uhbt, cs%vhbt, g, gv, us, &

939 cs%barotropic_CSp, cs%visc_rem_u, cs%visc_rem_v, spv_avg, cs%ADp, cs%OBC, cs%BT_cont, &

940 eta_pf_start, taux_bot, tauy_bot, uh_ptr, vh_ptr, u_ptr, v_ptr, etaav=eta_av)

941 if (cs%id_deta_dt>0) then

942 do j=js,je ; do i=is,ie ; deta_dt(i,j) = (eta_pred(i,j) - eta(i,j))*idt_bc ; enddo ; enddo

943 endif

944 do j=js,je ; do i=is,ie ; eta(i,j) = eta_pred(i,j) ; enddo ; enddo

945

946 call cpu_clock_end(id_clock_btstep)

947 if (showcalltree) call calltree_leave("btstep()")

948

949 if (cs%debug .and. debug_redundant) then

950 call check_redundant("u_accel_bt ", cs%u_accel_bt, cs%v_accel_bt, g, unscale=us%L_T2_to_m_s2)

951 endif

952

953 ! u = u + dt*( u_bc_accel + u_accel_bt )

954 call cpu_clock_begin(id_clock_mom_update)

955 !$OMP parallel do default(shared)

956 do k=1,nz

957 do j=js,je ; do i=isq,ieq

958 u_inst(i,j,k) = g%mask2dCu(i,j) * (u_inst(i,j,k) + dt * &

959 (u_bc_accel(i,j,k) + cs%u_accel_bt(i,j,k)))

960 enddo ; enddo

961 do j=jsq,jeq ; do i=is,ie

962 v_inst(i,j,k) = g%mask2dCv(i,j) * (v_inst(i,j,k) + dt * &

963 (v_bc_accel(i,j,k) + cs%v_accel_bt(i,j,k)))

964 enddo ; enddo

965 enddo

966 call cpu_clock_end(id_clock_mom_update)

967

968 if (cs%debug) then

969 call uvchksum("Corrector 1 [uv]", u_inst, v_inst, g%HI, haloshift=0, symmetric=sym, unscale=us%L_T_to_m_s)

970 call hchksum(h, "Corrector 1 h", g%HI, haloshift=1, unscale=gv%H_to_MKS)

971 call uvchksum("Corrector 1 [uv]h", uh, vh, g%HI, haloshift=2, &

972 symmetric=sym, unscale=gv%H_to_MKS*us%L_to_m**2*us%s_to_T)

973 ! call MOM_state_chksum("Corrector 1", u_inst, v_inst, h, uh, vh, G, GV, US, haloshift=1)

974 call mom_accel_chksum("Corrector accel", cs%CAu, cs%CAv, cs%PFu, cs%PFv, &

975 cs%diffu, cs%diffv, g, gv, us, cs%pbce, cs%u_accel_bt, cs%v_accel_bt, &

976 symmetric=sym)

977 endif

978

979 ! u <- u + dt d/dz visc d/dz u

980 ! u_av <- u_av + dt d/dz visc d/dz u_av

981 call cpu_clock_begin(id_clock_vertvisc)

982

983 ! if (CS%fpmix) then

984 ! uold(:,:,:) = 0.0

985 ! vold(:,:,:) = 0.0

986 ! do k=1,nz ; do j=js,je ; do I=Isq,Ieq

987 ! uold(I,j,k) = u_inst(I,j,k)

988 ! enddo ; enddo ; enddo

989 ! do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie

990 ! vold(i,J,k) = v_inst(i,J,k)

991 ! enddo ; enddo ; enddo

992 ! endif

993

994 call thickness_to_dz(h, tv, dz, g, gv, us, halo_size=1)

995 call vertvisc_coef(u_inst, v_inst, h, dz, forces, visc, tv, dt, g, gv, us, cs%vertvisc_CSp, cs%OBC, varmix)

996 call vertvisc(u_inst, v_inst, h, forces, visc, dt, cs%OBC, cs%ADp, cs%CDp, g, gv, us, &

997 cs%vertvisc_CSp, cs%taux_bot, cs%tauy_bot, waves=waves)

998

999 ! if (CS%fpmix) then

1000 ! call vertFPmix(u_inst, v_inst, uold, vold, hbl, h, forces, dt, &

1001 ! G, GV, US, CS%vertvisc_CSp, CS%OBC)

1002 ! call vertvisc(u_inst, v_inst, h, forces, visc, dt, CS%OBC, CS%ADp, CS%CDp, G, GV, US, &

1003 ! CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot, waves=waves)

1004 ! endif

1005

1006 if (g%nonblocking_updates) then

1007 call cpu_clock_end(id_clock_vertvisc)

1008 call start_group_pass(cs%pass_uv_inst, g%Domain, clock=id_clock_pass)

1009 call cpu_clock_begin(id_clock_vertvisc)

1010 endif

1011 call vertvisc_remnant(visc, cs%visc_rem_u, cs%visc_rem_v, dt, g, gv, us, cs%vertvisc_CSp)

1012 call cpu_clock_end(id_clock_vertvisc)

1013 if (showcalltree) call calltree_waypoint("done with vertvisc (step_MOM_dyn_split_RK2b)")

1014

1015 call do_group_pass(cs%pass_visc_rem, g%Domain, clock=id_clock_pass)

1016 if (g%nonblocking_updates) then

1017 call complete_group_pass(cs%pass_uv_inst, g%Domain, clock=id_clock_pass)

1018 else

1019 call do_group_pass(cs%pass_uv_inst, g%Domain, clock=id_clock_pass)

1020 endif

1021

1022 ! uh = u_av * h

1023 ! h = h + dt * div . uh

1024 ! u_av and v_av adjusted so their mass transports match uhbt and vhbt.

1025 call cpu_clock_begin(id_clock_continuity)

1026

1027 call continuity(u_inst, v_inst, h, h, uh, vh, dt, g, gv, us, cs%continuity_CSp, cs%OBC, pbv, &

1028 cs%uhbt, cs%vhbt, cs%visc_rem_u, cs%visc_rem_v, u_av, v_av, &

1029 du_cor=cs%du_av_inst, dv_cor=cs%dv_av_inst)

1030

1031 ! This tests the ability to readjust the instantaneous velocity, and here it changes answers only at roundoff.

1032 ! do k=1,nz ; do j=js,je ; do I=Isq,Ieq

1033 ! u_inst(I,j,k) = u_av(I,j,k) - CS%du_av_inst(I,j) * CS%visc_rem_u(I,j,k)

1034 ! enddo ; enddo ; enddo

1035 ! do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie

1036 ! v_inst(i,J,k) = v_av(i,J,k) - CS%dv_av_inst(i,J) * CS%visc_rem_v(i,J,k)

1037 ! enddo ; enddo ; enddo

1038

1039 call cpu_clock_end(id_clock_continuity)

1040

1041 call do_group_pass(cs%pass_h_uv, g%Domain, clock=id_clock_pass)

1042

1043 ! Whenever thickness changes let the diag manager know, target grids

1044 ! for vertical remapping may need to be regenerated.

1045 call diag_update_remap_grids(cs%diag)

1046 if (showcalltree) call calltree_waypoint("done with continuity (step_MOM_dyn_split_RK2b)")

1047

1048 if (associated(cs%OBC)) then

1049 call radiation_open_bdry_conds(cs%OBC, u_av, u_old_rad_obc, v_av, v_old_rad_obc, g, gv, us, dt)

1050 endif

1051

1052 !$OMP parallel do default(shared)

1053 do k=1,nz

1054 do j=js-2,je+2 ; do i=isq-2,ieq+2

1055 uhtr(i,j,k) = uhtr(i,j,k) + uh(i,j,k)*dt

1056 enddo ; enddo

1057 do j=jsq-2,jeq+2 ; do i=is-2,ie+2

1058 vhtr(i,j,k) = vhtr(i,j,k) + vh(i,j,k)*dt

1059 enddo ; enddo

1060 enddo

1061

1062 if (associated(cs%OBC)) then

1063 call update_segment_thickness_reservoirs(g, gv, uhtr, vhtr, h, cs%OBC)

1064 endif

1065

1066 ! if (CS%fpmix) then

1067 ! if (CS%id_uold > 0) call post_data(CS%id_uold, uold, CS%diag)

1068 ! if (CS%id_vold > 0) call post_data(CS%id_vold, vold, CS%diag)

1069 ! endif

1070

1071 ! The time-averaged free surface height has already been set by the last call to btstep.

1072

1073 ! Deallocate this memory to avoid a memory leak. ### We should revisit how this array is declared. -RWH

1074 if (dyn_p_surf .and. associated(eta_pf_start)) deallocate(eta_pf_start)

1075

1076 ! Here various terms used in to update the momentum equations are

1077 ! offered for time averaging.

1078 if (cs%id_PFu > 0) call post_data(cs%id_PFu, cs%PFu, cs%diag)

1079 if (cs%id_PFv > 0) call post_data(cs%id_PFv, cs%PFv, cs%diag)

1080 if (cs%id_CAu > 0) call post_data(cs%id_CAu, cs%CAu, cs%diag)

1081 if (cs%id_CAv > 0) call post_data(cs%id_CAv, cs%CAv, cs%diag)

1082

1083 ! Here the thickness fluxes are offered for time averaging.

1084 if (cs%id_uh > 0) call post_data(cs%id_uh, uh, cs%diag)

1085 if (cs%id_vh > 0) call post_data(cs%id_vh, vh, cs%diag)

1086 if (cs%id_uav > 0) call post_data(cs%id_uav, u_av, cs%diag)

1087 if (cs%id_vav > 0) call post_data(cs%id_vav, v_av, cs%diag)

1088 if (cs%id_u_BT_accel > 0) call post_data(cs%id_u_BT_accel, cs%u_accel_bt, cs%diag)

1089 if (cs%id_v_BT_accel > 0) call post_data(cs%id_v_BT_accel, cs%v_accel_bt, cs%diag)

1090

1091 ! Calculate effective areas and post data

1092 if (cs%id_ueffA > 0) then

1093 ueffa(:,:,:) = 0

1094 do k=1,nz ; do j=js,je ; do i=isq,ieq

1095 if (abs(u_av(i,j,k)) > 0.) ueffa(i,j,k) = uh(i,j,k) / u_av(i,j,k)

1096 enddo ; enddo ; enddo

1097 call post_data(cs%id_ueffA, ueffa, cs%diag)

1098 endif

1099

1100 if (cs%id_veffA > 0) then

1101 veffa(:,:,:) = 0

1102 do k=1,nz ; do j=jsq,jeq ; do i=is,ie

1103 if (abs(v_av(i,j,k)) > 0.) veffa(i,j,k) = vh(i,j,k) / v_av(i,j,k)

1104 enddo ; enddo ; enddo

1105 call post_data(cs%id_veffA, veffa, cs%diag)

1106 endif

1107

1108 ! Diagnostics of the fractional thicknesses times momentum budget terms

1109 ! 3D diagnostics hf_PFu etc. are commented because there is no clarity on proper remapping grid option.

1110 ! The code is retained for debugging purposes in the future.

1111 !if (CS%id_hf_PFu > 0) call post_product_u(CS%id_hf_PFu, CS%PFu, CS%ADp%diag_hfrac_u, G, nz, CS%diag)

1112 !if (CS%id_hf_PFv > 0) call post_product_v(CS%id_hf_PFv, CS%PFv, CS%ADp%diag_hfrac_v, G, nz, CS%diag)

1113 !if (CS%id_hf_CAu > 0) call post_product_u(CS%id_hf_CAu, CS%CAu, CS%ADp%diag_hfrac_u, G, nz, CS%diag)

1114 !if (CS%id_hf_CAv > 0) call post_product_v(CS%id_hf_CAv, CS%CAv, CS%ADp%diag_hfrac_v, G, nz, CS%diag)

1115 !if (CS%id_hf_u_BT_accel > 0) &

1116 ! call post_product_u(CS%id_hf_u_BT_accel, CS%u_accel_bt, CS%ADp%diag_hfrac_u, G, nz, CS%diag)

1117 !if (CS%id_hf_v_BT_accel > 0) &

1118 ! call post_product_v(CS%id_hf_v_BT_accel, CS%v_accel_bt, CS%ADp%diag_hfrac_v, G, nz, CS%diag)

1119

1120 ! Diagnostics for the vertical sum of layer thickness x prssure force accelerations

1121 if (cs%id_intz_PFu_2d > 0) call post_product_sum_u(cs%id_intz_PFu_2d, cs%PFu, cs%ADp%diag_hu, g, nz, cs%diag)

1122 if (cs%id_intz_PFv_2d > 0) call post_product_sum_v(cs%id_intz_PFv_2d, cs%PFv, cs%ADp%diag_hv, g, nz, cs%diag)

1123

1124 ! Diagnostics for thickness-weighted vertically averaged prssure force accelerations

1125 if (cs%id_hf_PFu_2d > 0) call post_product_sum_u(cs%id_hf_PFu_2d, cs%PFu, cs%ADp%diag_hfrac_u, g, nz, cs%diag)

1126 if (cs%id_hf_PFv_2d > 0) call post_product_sum_v(cs%id_hf_PFv_2d, cs%PFv, cs%ADp%diag_hfrac_v, g, nz, cs%diag)

1127

1128 ! Diagnostics for thickness x prssure force accelerations

1129 if (cs%id_h_PFu > 0) call post_product_u(cs%id_h_PFu, cs%PFu, cs%ADp%diag_hu, g, nz, cs%diag)

1130 if (cs%id_h_PFv > 0) call post_product_v(cs%id_h_PFv, cs%PFv, cs%ADp%diag_hv, g, nz, cs%diag)

1131

1132 ! Diagnostics of Coriolis acceleratations

1133 if (cs%id_intz_CAu_2d > 0) call post_product_sum_u(cs%id_intz_CAu_2d, cs%CAu, cs%ADp%diag_hu, g, nz, cs%diag)

1134 if (cs%id_intz_CAv_2d > 0) call post_product_sum_v(cs%id_intz_CAv_2d, cs%CAv, cs%ADp%diag_hv, g, nz, cs%diag)

1135 if (cs%id_hf_CAu_2d > 0) call post_product_sum_u(cs%id_hf_CAu_2d, cs%CAu, cs%ADp%diag_hfrac_u, g, nz, cs%diag)

1136 if (cs%id_hf_CAv_2d > 0) call post_product_sum_v(cs%id_hf_CAv_2d, cs%CAv, cs%ADp%diag_hfrac_v, g, nz, cs%diag)

1137 if (cs%id_h_CAu > 0) call post_product_u(cs%id_h_CAu, cs%CAu, cs%ADp%diag_hu, g, nz, cs%diag)

1138 if (cs%id_h_CAv > 0) call post_product_v(cs%id_h_CAv, cs%CAv, cs%ADp%diag_hv, g, nz, cs%diag)

1139

1140 ! Diagnostics of barotropic solver acceleratations

1141 if (cs%id_intz_u_BT_accel_2d > 0) &

1142 call post_product_sum_u(cs%id_intz_u_BT_accel_2d, cs%u_accel_bt, cs%ADp%diag_hu, g, nz, cs%diag)

1143 if (cs%id_intz_v_BT_accel_2d > 0) &

1144 call post_product_sum_v(cs%id_intz_v_BT_accel_2d, cs%v_accel_bt, cs%ADp%diag_hv, g, nz, cs%diag)

1145 if (cs%id_hf_u_BT_accel_2d > 0) &

1146 call post_product_sum_u(cs%id_hf_u_BT_accel_2d, cs%u_accel_bt, cs%ADp%diag_hfrac_u, g, nz, cs%diag)

1147 if (cs%id_hf_v_BT_accel_2d > 0) &

1148 call post_product_sum_v(cs%id_hf_v_BT_accel_2d, cs%v_accel_bt, cs%ADp%diag_hfrac_v, g, nz, cs%diag)

1149 if (cs%id_h_u_BT_accel > 0) &

1150 call post_product_u(cs%id_h_u_BT_accel, cs%u_accel_bt, cs%ADp%diag_hu, g, nz, cs%diag)

1151 if (cs%id_h_v_BT_accel > 0) &

1152 call post_product_v(cs%id_h_v_BT_accel, cs%v_accel_bt, cs%ADp%diag_hv, g, nz, cs%diag)

1153

1154 ! Diagnostics for momentum budget terms multiplied by visc_rem_[uv],

1155 if (cs%id_PFu_visc_rem > 0) call post_product_u(cs%id_PFu_visc_rem, cs%PFu, cs%ADp%visc_rem_u, g, nz, cs%diag)

1156 if (cs%id_PFv_visc_rem > 0) call post_product_v(cs%id_PFv_visc_rem, cs%PFv, cs%ADp%visc_rem_v, g, nz, cs%diag)

1157 if (cs%id_CAu_visc_rem > 0) call post_product_u(cs%id_CAu_visc_rem, cs%CAu, cs%ADp%visc_rem_u, g, nz, cs%diag)

1158 if (cs%id_CAv_visc_rem > 0) call post_product_v(cs%id_CAv_visc_rem, cs%CAv, cs%ADp%visc_rem_v, g, nz, cs%diag)

1159 if (cs%id_u_BT_accel_visc_rem > 0) &

1160 call post_product_u(cs%id_u_BT_accel_visc_rem, cs%u_accel_bt, cs%ADp%visc_rem_u, g, nz, cs%diag)

1161 if (cs%id_v_BT_accel_visc_rem > 0) &

1162 call post_product_v(cs%id_v_BT_accel_visc_rem, cs%v_accel_bt, cs%ADp%visc_rem_v, g, nz, cs%diag)

1163

1164 ! Diagnostics related to changes in eta

1165 if (cs%id_deta_dt > 0) call post_data(cs%id_deta_dt, deta_dt, cs%diag)

1166

1167 if (cs%debug) then

1168 call mom_state_chksum("Corrector ", u_av, v_av, h, uh, vh, g, gv, us, symmetric=sym)

1169 ! call uvchksum("Corrector inst [uv]", u_inst, v_inst, G%HI, symmetric=sym, unscale=US%L_T_to_m_s)

1170 endif

1171

1172 if (showcalltree) call calltree_leave("step_MOM_dyn_split_RK2b()")

1173

1174end subroutine step_mom_dyn_split_rk2b

1175

1176!> This subroutine sets up any auxiliary restart variables that are specific

1177!! to the split-explicit time stepping scheme. All variables registered here should

1178!! have the ability to be recreated if they are not present in a restart file.

1179subroutine register_restarts_dyn_split_rk2b(HI, GV, US, param_file, CS, restart_CS, uh, vh)

1180 type(hor_index_type), intent(in) :: hi !< Horizontal index structure

1181 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure

1182 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type

1183 type(param_file_type), intent(in) :: param_file !< parameter file

1184 type(mom_dyn_split_rk2b_cs), pointer :: cs !< module control structure

1185 type(mom_restart_cs), intent(inout) :: restart_cs !< MOM restart control structure

1186 real, dimension(SZIB_(HI),SZJ_(HI),SZK_(GV)), &

1187 target, intent(inout) :: uh !< zonal volume or mass transport [H L2 T-1 ~> m3 s-1 or kg s-1]

1188 real, dimension(SZI_(HI),SZJB_(HI),SZK_(GV)), &

1189 target, intent(inout) :: vh !< merid volume or mass transport [H L2 T-1 ~> m3 s-1 or kg s-1]

1190

1191 type(vardesc) :: vd(2)

1192 character(len=48) :: thickness_units, flux_units

1193

1194 integer :: isd, ied, jsd, jed, nz, isdb, iedb, jsdb, jedb

1195

1196 isd = hi%isd ; ied = hi%ied ; jsd = hi%jsd ; jed = hi%jed ; nz = gv%ke

1197 isdb = hi%IsdB ; iedb = hi%IedB ; jsdb = hi%JsdB ; jedb = hi%JedB

1198

1199 ! This is where a control structure specific to this module would be allocated.

1200 if (associated(cs)) then

1201 call mom_error(warning, "register_restarts_dyn_split_RK2b called with an associated "// &

1202 "control structure.")

1203 return

1204 endif

1205 allocate(cs)

1206

1207 alloc_(cs%diffu(isdb:iedb,jsd:jed,nz)) ; cs%diffu(:,:,:) = 0.0

1208 alloc_(cs%diffv(isd:ied,jsdb:jedb,nz)) ; cs%diffv(:,:,:) = 0.0

1209 alloc_(cs%CAu(isdb:iedb,jsd:jed,nz)) ; cs%CAu(:,:,:) = 0.0

1210 alloc_(cs%CAv(isd:ied,jsdb:jedb,nz)) ; cs%CAv(:,:,:) = 0.0

1211 alloc_(cs%CAu_pred(isdb:iedb,jsd:jed,nz)) ; cs%CAu_pred(:,:,:) = 0.0

1212 alloc_(cs%CAv_pred(isd:ied,jsdb:jedb,nz)) ; cs%CAv_pred(:,:,:) = 0.0

1213 alloc_(cs%PFu(isdb:iedb,jsd:jed,nz)) ; cs%PFu(:,:,:) = 0.0

1214 alloc_(cs%PFv(isd:ied,jsdb:jedb,nz)) ; cs%PFv(:,:,:) = 0.0

1215 alloc_(cs%du_av_inst(isdb:iedb,jsd:jed)) ; cs%du_av_inst(:,:) = 0.0

1216 alloc_(cs%dv_av_inst(isd:ied,jsdb:jedb)) ; cs%dv_av_inst(:,:) = 0.0

1217

1218 alloc_(cs%eta(isd:ied,jsd:jed)) ; cs%eta(:,:) = 0.0

1219

1220 thickness_units = get_thickness_units(gv)

1221 flux_units = get_flux_units(gv)

1222

1223 if (gv%Boussinesq) then

1224 call register_restart_field(cs%eta, "sfc", .false., restart_cs, &

1225 longname="Free surface Height", units=thickness_units, conversion=gv%H_to_mks)

1226 else

1227 call register_restart_field(cs%eta, "p_bot", .false., restart_cs, &

1228 longname="Bottom Pressure", units=thickness_units, conversion=gv%H_to_mks)

1229 endif

1230

1231 ! These are needed to reconstruct the phase in the barotorpic solution.

1232 vd(1) = var_desc("du_avg_inst", "m s-1", &

1233 "Barotropic velocity increment between instantaneous and filtered zonal velocities", 'u', '1')

1234 vd(2) = var_desc("dv_avg_inst", "m s-1", &

1235 "Barotropic velocity increment between instantaneous and filtered meridional velocities", 'v', '1')

1236 call register_restart_pair(cs%du_av_inst, cs%dv_av_inst, vd(1), vd(2), .false., restart_cs, &

1237 conversion=us%L_T_to_m_s)

1238

1239 call register_barotropic_restarts(hi, gv, us, param_file, cs%barotropic_CSp, restart_cs)

1240

1241end subroutine register_restarts_dyn_split_rk2b

1242

1243!> This subroutine does remapping for the auxiliary restart variables that are used

1244!! with the split RK2 time stepping scheme.

1245subroutine remap_dyn_split_rk2b_aux_vars(G, GV, CS, h_old_u, h_old_v, h_new_u, h_new_v, ALE_CSp)

1246 type(ocean_grid_type), intent(inout) :: g !< ocean grid structure

1247 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure

1248 type(mom_dyn_split_rk2b_cs), pointer :: cs !< module control structure

1249 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &

1250 intent(in) :: h_old_u !< Source grid thickness at zonal

1251 !! velocity points [H ~> m or kg m-2]

1252 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &

1253 intent(in) :: h_old_v !< Source grid thickness at meridional

1254 !! velocity points [H ~> m or kg m-2]

1255 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &

1256 intent(in) :: h_new_u !< Destination grid thickness at zonal

1257 !! velocity points [H ~> m or kg m-2]

1258 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &

1259 intent(in) :: h_new_v !< Destination grid thickness at meridional

1260 !! velocity points [H ~> m or kg m-2]

1261 type(ale_cs), pointer :: ale_csp !< ALE control structure to use when remapping

1262

1263 return

1264

1265end subroutine remap_dyn_split_rk2b_aux_vars

1266

1267!> This subroutine initializes all of the variables that are used by this

1268!! dynamic core, including diagnostics and the cpu clocks.

1269subroutine initialize_dyn_split_rk2b(u, v, h, tv, uh, vh, eta, Time, G, GV, US, param_file, &

1270 diag, CS, HA_CSp, restart_CS, dt, Accel_diag, Cont_diag, MIS, &

1271 VarMix, MEKE, thickness_diffuse_CSp, &

1272 OBC, update_OBC_CSp, ALE_CSp, set_visc, &

1273 visc, dirs, ntrunc, pbv, calc_dtbt, cont_stencil, dyn_h_stencil)

1274 type(ocean_grid_type), intent(inout) :: g !< ocean grid structure

1275 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure

1276 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type

1277 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &

1278 intent(inout) :: u !< zonal velocity [L T-1 ~> m s-1]

1279 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &

1280 intent(inout) :: v !< merid velocity [L T-1 ~> m s-1]

1281 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

1282 intent(inout) :: h !< layer thickness [H ~> m or kg m-2]

1283 type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic type

1284 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &

1285 target, intent(inout) :: uh !< zonal volume/mass transport [H L2 T-1 ~> m3 s-1 or kg s-1]

1286 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &

1287 target, intent(inout) :: vh !< merid volume/mass transport [H L2 T-1 ~> m3 s-1 or kg s-1]

1288 real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: eta !< free surface height or column mass [H ~> m or kg m-2]

1289 type(time_type), target, intent(in) :: time !< current model time

1290 type(param_file_type), intent(in) :: param_file !< parameter file for parsing

1291 type(diag_ctrl), target, intent(inout) :: diag !< to control diagnostics

1292 type(mom_dyn_split_rk2b_cs), pointer :: cs !< module control structure

1293 type(harmonic_analysis_cs), pointer :: ha_csp !< A pointer to the control structure of the

1294 !! harmonic analysis module

1295 type(mom_restart_cs), intent(inout) :: restart_cs !< MOM restart control structure

1296 real, intent(in) :: dt !< time step [T ~> s]

1297 type(accel_diag_ptrs), target, intent(inout) :: accel_diag !< points to momentum equation terms for

1298 !! budget analysis

1299 type(cont_diag_ptrs), target, intent(inout) :: cont_diag !< points to terms in continuity equation

1300 type(ocean_internal_state), intent(inout) :: mis !< "MOM6 internal state" used to pass

1301 !! diagnostic pointers

1302 type(varmix_cs), intent(inout) :: varmix !< points to spatially variable viscosities

1303 type(meke_type), intent(inout) :: meke !< MEKE fields

1304 type(thickness_diffuse_cs), intent(inout) :: thickness_diffuse_csp !< Pointer to the control structure

1305 !! used for the isopycnal height diffusive transport.

1306 type(ocean_obc_type), pointer :: obc !< points to OBC related fields

1307 type(update_obc_cs), pointer :: update_obc_csp !< points to OBC update related fields

1308 type(ale_cs), pointer :: ale_csp !< points to ALE control structure

1309 type(set_visc_cs), target, intent(in) :: set_visc !< set_visc control structure

1310 type(vertvisc_type), intent(inout) :: visc !< vertical viscosities, bottom drag, and related

1311 type(directories), intent(in) :: dirs !< contains directory paths

1312 integer, target, intent(inout) :: ntrunc !< A target for the variable that records

1313 !! the number of times the velocity is

1314 !! truncated (this should be 0).

1315 logical, intent(out) :: calc_dtbt !< If true, recalculate the barotropic time step

1316 type(porous_barrier_type), intent(in) :: pbv !< porous barrier fractional cell metrics

1317 integer, intent(out) :: cont_stencil !< The stencil for thickness

1318 !! from the continuity solver.

1319 integer, intent(out) :: dyn_h_stencil !< The stencil for thickness for the

1320 !! dynamics based on the continuity

1321 !! solver and Coriolis scheme.

1322

1323 ! local variables

1324 character(len=40) :: mdl = "MOM_dynamics_split_RK2b" ! This module's name.

1325 ! This include declares and sets the variable "version".

1326# include "version_variable.h"

1327 character(len=48) :: thickness_units, flux_units, eta_rest_name

1328 logical :: debug_truncations

1329 logical :: enable_bugs ! If true, the defaults for recently added bug-fix flags are set to

1330 ! recreate the bugs, or if false bugs are only used if actively selected.

1331 logical :: visc_rem_bug ! Stores the value of runtime paramter VISC_REM_BUG.

1332

1333 integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz

1334 integer :: isdb, iedb, jsdb, jedb

1335 integer :: nc ! Number of tidal constituents to be harmonically analyzed

1336 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke

1337 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed

1338 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB

1339

1340 if (.not.associated(cs)) call mom_error(fatal, &

1341 "initialize_dyn_split_RK2b called with an unassociated control structure.")

1342 if (cs%module_is_initialized) then

1343 call mom_error(warning, "initialize_dyn_split_RK2b called with a control "// &

1344 "structure that has already been initialized.")

1345 return

1346 endif

1347 cs%module_is_initialized = .true.

1348

1349 cs%diag => diag

1350

1351 call log_version(param_file, mdl, version, "")

1352 call get_param(param_file, mdl, "TIDES", cs%use_tides, &

1353 "If true, apply tidal momentum forcing.", default=.false.)

1354 call get_param(param_file, mdl, "CALCULATE_SAL", cs%calculate_SAL, &

1355 "If true, calculate self-attraction and loading.", default=cs%use_tides)

1356 call get_param(param_file, mdl, "USE_HA", cs%use_HA, &

1357 "If true, perform inline harmonic analysis.", default=.false.)

1358 call get_param(param_file, mdl, "HA_N_CONST", nc, &

1359 "Number of tidal constituents to be harmonically analyzed.", &

1360 default=0, do_not_log=.not.cs%use_HA)

1361 if (nc<=0) cs%use_HA = .false.

1362 call get_param(param_file, mdl, "BE", cs%be, &

1363 "If SPLIT is true, BE determines the relative weighting "//&

1364 "of a 2nd-order Runga-Kutta baroclinic time stepping "//&

1365 "scheme (0.5) and a backward Euler scheme (1) that is "//&

1366 "used for the Coriolis and inertial terms. BE may be "//&

1367 "from 0.5 to 1, but instability may occur near 0.5. "//&

1368 "BE is also applicable if SPLIT is false and USE_RK2 "//&

1369 "is true.", units="nondim", default=0.6)

1370 call get_param(param_file, mdl, "BEGW", cs%begw, &

1371 "If SPLIT is true, BEGW is a number from 0 to 1 that "//&

1372 "controls the extent to which the treatment of gravity "//&

1373 "waves is forward-backward (0) or simulated backward "//&

1374 "Euler (1). 0 is almost always used. "//&

1375 "If SPLIT is false and USE_RK2 is true, BEGW can be "//&

1376 "between 0 and 0.5 to damp gravity waves.", &

1377 units="nondim", default=0.0)

1378 call get_param(param_file, mdl, "SET_DTBT_USE_BT_CONT", cs%dtbt_use_bt_cont, &

1379 "If true, use BT_CONT to calculate DTBT if possible.", default=.false.)

1380 call get_param(param_file, mdl, "SPLIT_BOTTOM_STRESS", cs%split_bottom_stress, &

1381 "If true, provide the bottom stress calculated by the "//&

1382 "vertical viscosity to the barotropic solver.", default=.false.)

1383 call get_param(param_file, mdl, "BT_ADJ_CORR_MASS_SRC", cs%BT_adj_corr_mass_src, &

1384 "If true, recalculates the barotropic mass source after "//&

1385 "predictor step. This should make little difference in the "//&

1386 "deep ocean but appears to help for vanished layers. If false, "//&

1387 "uses the same mass source as from the predictor step.", default=.true.)

1388 ! call get_param(param_file, mdl, "FPMIX", CS%fpmix, &

1389 ! "If true, apply profiles of momentum flux magnitude and direction.", &

1390 ! default=.false.)

1391 cs%fpmix = .false.

1392 call get_param(param_file, mdl, "REMAP_AUXILIARY_VARS", cs%remap_aux, &

1393 "If true, apply ALE remapping to all of the auxiliary 3-dimensional "//&

1394 "variables that are needed to reproduce across restarts, similarly to "//&

1395 "what is already being done with the primary state variables. "//&

1396 "The default should be changed to true.", default=.false., do_not_log=.true.)

1397 call get_param(param_file, mdl, "DEBUG", cs%debug, &

1398 "If true, write out verbose debugging data.", &

1399 default=.false., debuggingparam=.true.)

1400 call get_param(param_file, mdl, "OBC_DEBUGGING_TESTS", cs%debug_OBC, &

1401 "If true, do additional calls resetting certain values to help verify the "//&

1402 "correctness of the open boundary condition code.", &

1403 default=.false., old_name="DEBUG_OBC", debuggingparam=.true., do_not_log=.true.)

1404 call get_param(param_file, mdl, "DEBUG_TRUNCATIONS", debug_truncations, &

1405 default=.false.)

1406 call get_param(param_file, mdl, "ENABLE_BUGS_BY_DEFAULT", enable_bugs, &

1407 default=.true., do_not_log=.true.) ! This is logged from MOM.F90.

1408 call get_param(param_file, mdl, "VISC_REM_BUG", visc_rem_bug, &

1409 "If true, visc_rem_[uv] in split mode is incorrectly calculated or accounted "//&

1410 "for in two places. This parameter controls the defaults of two individual "//&

1411 "flags, VISC_REM_TIMESTEP_BUG in MOM_dynamics_split_RK2(b) and "//&

1412 "VISC_REM_BT_WEIGHT_BUG in MOM_barotropic.", default=.false.)

1413 call get_param(param_file, mdl, "VISC_REM_TIMESTEP_BUG", cs%visc_rem_dt_bug, &

1414 "If true, recover a bug that uses dt_pred rather than dt in "//&

1415 "vertvisc_remnant() at the end of predictor stage for the following "//&

1416 "continuity() and btstep() calls in the corrector step. Default of this flag "//&

1417 "is set by VISC_REM_BUG", default=visc_rem_bug)

1418

1419

1420 allocate(cs%taux_bot(isdb:iedb,jsd:jed), source=0.0)

1421 allocate(cs%tauy_bot(isd:ied,jsdb:jedb), source=0.0)

1422

1423 alloc_(cs%uhbt(isdb:iedb,jsd:jed)) ; cs%uhbt(:,:) = 0.0

1424 alloc_(cs%vhbt(isd:ied,jsdb:jedb)) ; cs%vhbt(:,:) = 0.0

1425 alloc_(cs%visc_rem_u(isdb:iedb,jsd:jed,nz)) ; cs%visc_rem_u(:,:,:) = 0.0

1426 alloc_(cs%visc_rem_v(isd:ied,jsdb:jedb,nz)) ; cs%visc_rem_v(:,:,:) = 0.0

1427 alloc_(cs%eta_PF(isd:ied,jsd:jed)) ; cs%eta_PF(:,:) = 0.0

1428 alloc_(cs%pbce(isd:ied,jsd:jed,nz)) ; cs%pbce(:,:,:) = 0.0

1429

1430 alloc_(cs%u_accel_bt(isdb:iedb,jsd:jed,nz)) ; cs%u_accel_bt(:,:,:) = 0.0

1431 alloc_(cs%v_accel_bt(isd:ied,jsdb:jedb,nz)) ; cs%v_accel_bt(:,:,:) = 0.0

1432 alloc_(cs%PFu_Stokes(isdb:iedb,jsd:jed,nz)) ; cs%PFu_Stokes(:,:,:) = 0.0

1433 alloc_(cs%PFv_Stokes(isd:ied,jsdb:jedb,nz)) ; cs%PFv_Stokes(:,:,:) = 0.0

1434

1435 mis%diffu => cs%diffu

1436 mis%diffv => cs%diffv

1437 mis%PFu => cs%PFu

1438 mis%PFv => cs%PFv

1439 mis%CAu => cs%CAu

1440 mis%CAv => cs%CAv

1441 mis%pbce => cs%pbce

1442 mis%u_accel_bt => cs%u_accel_bt

1443 mis%v_accel_bt => cs%v_accel_bt

1444

1445 cs%ADp => accel_diag

1446 cs%CDp => cont_diag

1447 accel_diag%diffu => cs%diffu

1448 accel_diag%diffv => cs%diffv

1449 accel_diag%PFu => cs%PFu

1450 accel_diag%PFv => cs%PFv

1451 accel_diag%CAu => cs%CAu

1452 accel_diag%CAv => cs%CAv

1453 accel_diag%u_accel_bt => cs%u_accel_bt

1454 accel_diag%v_accel_bt => cs%v_accel_bt

1455

1456 allocate(cs%AD_pred)

1457 cs%AD_pred%diffu => cs%diffu

1458 cs%AD_pred%diffv => cs%diffv

1459 cs%AD_pred%PFu => cs%PFu

1460 cs%AD_pred%PFv => cs%PFv

1461 cs%AD_pred%CAu => cs%CAu_pred

1462 cs%AD_pred%CAv => cs%CAv_pred

1463 cs%AD_pred%u_accel_bt => cs%u_accel_bt

1464 cs%AD_pred%v_accel_bt => cs%v_accel_bt

1465

1466! Accel_diag%pbce => CS%pbce

1467! Accel_diag%u_accel_bt => CS%u_accel_bt ; Accel_diag%v_accel_bt => CS%v_accel_bt

1468! Accel_diag%u_av => CS%u_av ; Accel_diag%v_av => CS%v_av

1469

1470 call continuity_init(time, g, gv, us, param_file, diag, cs%continuity_CSp, cs%OBC)

1471 cont_stencil = continuity_stencil(cs%continuity_CSp)

1472 call coriolisadv_init(time, g, gv, us, param_file, diag, cs%ADp, cs%CoriolisAdv)

1473 dyn_h_stencil = max(cont_stencil, coriolisadv_stencil(cs%CoriolisAdv))

1474 if (cs%calculate_SAL) call sal_init(h, tv, g, gv, us, param_file, cs%SAL_CSp, restart_cs)

1475 if (cs%use_tides) call tidal_forcing_init(time, g, us, param_file, cs%tides_CSp)

1476 if (cs%use_HA) then

1477 call ha_init(time, us, param_file, nc, cs%HA_CSp)

1478 ha_csp => cs%HA_CSp

1479 else

1480 ha_csp => null()

1481 endif

1482 call pressureforce_init(time, g, gv, us, param_file, diag, cs%PressureForce_CSp, cs%ADp, &

1483 cs%SAL_CSp, cs%tides_CSp)

1484 call hor_visc_init(time, g, gv, us, param_file, diag, cs%hor_visc, adp=cs%ADp)

1485 call vertvisc_init(mis, time, g, gv, us, param_file, diag, cs%ADp, dirs, &

1486 ntrunc, cs%vertvisc_CSp)

1487 cs%set_visc_CSp => set_visc

1488 call updatecfltruncationvalue(time, cs%vertvisc_CSp, us, &

1489 activate=is_new_run(restart_cs) )

1490

1491 if (associated(ale_csp)) cs%ALE_CSp => ale_csp

1492 if (associated(obc)) then

1493 cs%OBC => obc

1494 if (obc%ramp) call update_obc_ramp(time, cs%OBC, us, &

1495 activate=is_new_run(restart_cs) )

1496 endif

1497 if (associated(update_obc_csp)) cs%update_OBC_CSp => update_obc_csp

1498

1499 eta_rest_name = "sfc" ; if (.not.gv%Boussinesq) eta_rest_name = "p_bot"

1500 if (.not. query_initialized(cs%eta, trim(eta_rest_name), restart_cs)) then

1501 ! Estimate eta based on the layer thicknesses - h. With the Boussinesq

1502 ! approximation, eta is the free surface height anomaly, while without it

1503 ! eta is the mass of ocean per unit area. eta always has the same

1504 ! dimensions as h, either m or kg m-3.

1505 ! CS%eta(:,:) = 0.0 already from initialization.

1506 if (gv%Boussinesq) then

1507 do j=js,je ; do i=is,ie ; cs%eta(i,j) = -gv%Z_to_H * g%bathyT(i,j) ; enddo ; enddo

1508 endif

1509 do k=1,nz ; do j=js,je ; do i=is,ie

1510 cs%eta(i,j) = cs%eta(i,j) + h(i,j,k)

1511 enddo ; enddo ; enddo

1512 call set_initialized(cs%eta, trim(eta_rest_name), restart_cs)

1513 endif

1514 ! Copy eta into an output array.

1515 do j=js,je ; do i=is,ie ; eta(i,j) = cs%eta(i,j) ; enddo ; enddo

1516

1517 call barotropic_init(u, v, h, time, g, gv, us, param_file, diag, &

1518 cs%barotropic_CSp, restart_cs, calc_dtbt, cs%BT_cont, &

1519 cs%OBC, cs%SAL_CSp, ha_csp)

1520

1521 flux_units = get_flux_units(gv)

1522 thickness_units = get_thickness_units(gv)

1523 cs%id_uh = register_diag_field('ocean_model', 'uh', diag%axesCuL, time, &

1524 'Zonal Thickness Flux', flux_units, conversion=gv%H_to_MKS*us%L_to_m**2*us%s_to_T, &

1525 y_cell_method='sum', v_extensive=.true.)

1526 cs%id_vh = register_diag_field('ocean_model', 'vh', diag%axesCvL, time, &

1527 'Meridional Thickness Flux', flux_units, conversion=gv%H_to_MKS*us%L_to_m**2*us%s_to_T, &

1528 x_cell_method='sum', v_extensive=.true.)

1529

1530 cs%id_CAu = register_diag_field('ocean_model', 'CAu', diag%axesCuL, time, &

1531 'Zonal Coriolis and Advective Acceleration', 'm s-2', conversion=us%L_T2_to_m_s2)

1532 cs%id_CAv = register_diag_field('ocean_model', 'CAv', diag%axesCvL, time, &

1533 'Meridional Coriolis and Advective Acceleration', 'm s-2', conversion=us%L_T2_to_m_s2)

1534 cs%id_PFu = register_diag_field('ocean_model', 'PFu', diag%axesCuL, time, &

1535 'Zonal Pressure Force Acceleration', 'm s-2', conversion=us%L_T2_to_m_s2)

1536 cs%id_PFv = register_diag_field('ocean_model', 'PFv', diag%axesCvL, time, &

1537 'Meridional Pressure Force Acceleration', 'm s-2', conversion=us%L_T2_to_m_s2)

1538 cs%id_ueffA = register_diag_field('ocean_model', 'ueffA', diag%axesCuL, time, &

1539 'Effective U-Face Area', 'm^2', conversion=gv%H_to_m*us%L_to_m, &

1540 y_cell_method='sum', v_extensive=.true.)

1541 cs%id_veffA = register_diag_field('ocean_model', 'veffA', diag%axesCvL, time, &

1542 'Effective V-Face Area', 'm^2', conversion=gv%H_to_m*us%L_to_m, &

1543 x_cell_method='sum', v_extensive=.true.)

1544 if (gv%Boussinesq) then

1545 cs%id_deta_dt = register_diag_field('ocean_model', 'deta_dt', diag%axesT1, time, &

1546 'Barotropic SSH tendency due to dynamics', trim(thickness_units)//' s-1', conversion=gv%H_to_MKS*us%s_to_T)

1547 else

1548 cs%id_deta_dt = register_diag_field('ocean_model', 'deta_dt', diag%axesT1, time, &

1549 'Barotropic column-mass tendency due to dynamics', trim(thickness_units)//' s-1', &

1550 conversion=gv%H_to_mks*us%s_to_T)

1551 endif

1552

1553 !CS%id_hf_PFu = register_diag_field('ocean_model', 'hf_PFu', diag%axesCuL, Time, &

1554 ! 'Fractional Thickness-weighted Zonal Pressure Force Acceleration', &

1555 ! 'm s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)

1556 !if (CS%id_hf_PFu > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_u,IsdB,IedB,jsd,jed,nz)

1557

1558 !CS%id_hf_PFv = register_diag_field('ocean_model', 'hf_PFv', diag%axesCvL, Time, &

1559 ! 'Fractional Thickness-weighted Meridional Pressure Force Acceleration', &

1560 ! 'm s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)

1561 !if (CS%id_hf_PFv > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz)

1562

1563 !CS%id_hf_CAu = register_diag_field('ocean_model', 'hf_CAu', diag%axesCuL, Time, &

1564 ! 'Fractional Thickness-weighted Zonal Coriolis and Advective Acceleration', &

1565 ! 'm s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)

1566 !if (CS%id_hf_CAu > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_u,IsdB,IedB,jsd,jed,nz)

1567

1568 !CS%id_hf_CAv = register_diag_field('ocean_model', 'hf_CAv', diag%axesCvL, Time, &

1569 ! 'Fractional Thickness-weighted Meridional Coriolis and Advective Acceleration', &

1570 ! 'm s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)

1571 !if (CS%id_hf_CAv > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz)

1572

1573 cs%id_hf_PFu_2d = register_diag_field('ocean_model', 'hf_PFu_2d', diag%axesCu1, time, &

1574 'Depth-sum Fractional Thickness-weighted Zonal Pressure Force Acceleration', &

1575 'm s-2', conversion=us%L_T2_to_m_s2)

1576 if (cs%id_hf_PFu_2d > 0) call safe_alloc_ptr(cs%ADp%diag_hfrac_u,isdb,iedb,jsd,jed,nz)

1577

1578 cs%id_hf_PFv_2d = register_diag_field('ocean_model', 'hf_PFv_2d', diag%axesCv1, time, &

1579 'Depth-sum Fractional Thickness-weighted Meridional Pressure Force Acceleration', &

1580 'm s-2', conversion=us%L_T2_to_m_s2)

1581 if (cs%id_hf_PFv_2d > 0) call safe_alloc_ptr(cs%ADp%diag_hfrac_v,isd,ied,jsdb,jedb,nz)

1582

1583 cs%id_h_PFu = register_diag_field('ocean_model', 'h_PFu', diag%axesCuL, time, &

1584 'Thickness Multiplied Zonal Pressure Force Acceleration', &

1585 'm2 s-2', conversion=gv%H_to_m*us%L_T2_to_m_s2)

1586 if (cs%id_h_PFu > 0) call safe_alloc_ptr(cs%ADp%diag_hu,isdb,iedb,jsd,jed,nz)

1587

1588 cs%id_h_PFv = register_diag_field('ocean_model', 'h_PFv', diag%axesCvL, time, &

1589 'Thickness Multiplied Meridional Pressure Force Acceleration', &

1590 'm2 s-2', conversion=gv%H_to_m*us%L_T2_to_m_s2)

1591 if (cs%id_h_PFv > 0) call safe_alloc_ptr(cs%ADp%diag_hv,isd,ied,jsdb,jedb,nz)

1592

1593 cs%id_intz_PFu_2d = register_diag_field('ocean_model', 'intz_PFu_2d', diag%axesCu1, time, &

1594 'Depth-integral of Zonal Pressure Force Acceleration', &

1595 'm2 s-2', conversion=gv%H_to_m*us%L_T2_to_m_s2)

1596 if (cs%id_intz_PFu_2d > 0) call safe_alloc_ptr(cs%ADp%diag_hu,isdb,iedb,jsd,jed,nz)

1597

1598 cs%id_intz_PFv_2d = register_diag_field('ocean_model', 'intz_PFv_2d', diag%axesCv1, time, &

1599 'Depth-integral of Meridional Pressure Force Acceleration', &

1600 'm2 s-2', conversion=gv%H_to_m*us%L_T2_to_m_s2)

1601 if (cs%id_intz_PFv_2d > 0) call safe_alloc_ptr(cs%ADp%diag_hv,isd,ied,jsdb,jedb,nz)

1602

1603 cs%id_hf_CAu_2d = register_diag_field('ocean_model', 'hf_CAu_2d', diag%axesCu1, time, &

1604 'Depth-sum Fractional Thickness-weighted Zonal Coriolis and Advective Acceleration', &

1605 'm s-2', conversion=us%L_T2_to_m_s2)

1606 if (cs%id_hf_CAu_2d > 0) call safe_alloc_ptr(cs%ADp%diag_hfrac_u,isdb,iedb,jsd,jed,nz)

1607

1608 cs%id_hf_CAv_2d = register_diag_field('ocean_model', 'hf_CAv_2d', diag%axesCv1, time, &

1609 'Depth-sum Fractional Thickness-weighted Meridional Coriolis and Advective Acceleration', &

1610 'm s-2', conversion=us%L_T2_to_m_s2)

1611 if (cs%id_hf_CAv_2d > 0) call safe_alloc_ptr(cs%ADp%diag_hfrac_v,isd,ied,jsdb,jedb,nz)

1612

1613 cs%id_h_CAu = register_diag_field('ocean_model', 'h_CAu', diag%axesCuL, time, &

1614 'Thickness Multiplied Zonal Coriolis and Advective Acceleration', &

1615 'm2 s-2', conversion=gv%H_to_m*us%L_T2_to_m_s2)

1616 if (cs%id_h_CAu > 0) call safe_alloc_ptr(cs%ADp%diag_hu,isdb,iedb,jsd,jed,nz)

1617

1618 cs%id_h_CAv = register_diag_field('ocean_model', 'h_CAv', diag%axesCvL, time, &

1619 'Thickness Multiplied Meridional Coriolis and Advective Acceleration', &

1620 'm2 s-2', conversion=gv%H_to_m*us%L_T2_to_m_s2)

1621 if (cs%id_h_CAv > 0) call safe_alloc_ptr(cs%ADp%diag_hv,isd,ied,jsdb,jedb,nz)

1622

1623 cs%id_intz_CAu_2d = register_diag_field('ocean_model', 'intz_CAu_2d', diag%axesCu1, time, &

1624 'Depth-integral of Zonal Coriolis and Advective Acceleration', &

1625 'm2 s-2', conversion=gv%H_to_m*us%L_T2_to_m_s2)

1626 if (cs%id_intz_CAu_2d > 0) call safe_alloc_ptr(cs%ADp%diag_hu,isdb,iedb,jsd,jed,nz)

1627

1628 cs%id_intz_CAv_2d = register_diag_field('ocean_model', 'intz_CAv_2d', diag%axesCv1, time, &

1629 'Depth-integral of Meridional Coriolis and Advective Acceleration', &

1630 'm2 s-2', conversion=gv%H_to_m*us%L_T2_to_m_s2)

1631 if (cs%id_intz_CAv_2d > 0) call safe_alloc_ptr(cs%ADp%diag_hv,isd,ied,jsdb,jedb,nz)

1632

1633 cs%id_uav = register_diag_field('ocean_model', 'uav', diag%axesCuL, time, &

1634 'Barotropic-step Averaged Zonal Velocity', 'm s-1', conversion=us%L_T_to_m_s)

1635 cs%id_vav = register_diag_field('ocean_model', 'vav', diag%axesCvL, time, &

1636 'Barotropic-step Averaged Meridional Velocity', 'm s-1', conversion=us%L_T_to_m_s)

1637

1638 cs%id_u_BT_accel = register_diag_field('ocean_model', 'u_BT_accel', diag%axesCuL, time, &

1639 'Barotropic Anomaly Zonal Acceleration', 'm s-2', conversion=us%L_T2_to_m_s2)

1640 cs%id_v_BT_accel = register_diag_field('ocean_model', 'v_BT_accel', diag%axesCvL, time, &

1641 'Barotropic Anomaly Meridional Acceleration', 'm s-2', conversion=us%L_T2_to_m_s2)

1642

1643 !CS%id_hf_u_BT_accel = register_diag_field('ocean_model', 'hf_u_BT_accel', diag%axesCuL, Time, &

1644 ! 'Fractional Thickness-weighted Barotropic Anomaly Zonal Acceleration', &

1645 ! 'm s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)

1646 !if (CS%id_hf_u_BT_accel > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_u,IsdB,IedB,jsd,jed,nz)

1647

1648 !CS%id_hf_v_BT_accel = register_diag_field('ocean_model', 'hf_v_BT_accel', diag%axesCvL, Time, &

1649 ! 'Fractional Thickness-weighted Barotropic Anomaly Meridional Acceleration', &

1650 ! 'm s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)

1651 !if (CS%id_hf_v_BT_accel > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz)

1652

1653 cs%id_hf_u_BT_accel_2d = register_diag_field('ocean_model', 'hf_u_BT_accel_2d', diag%axesCu1, time, &

1654 'Depth-sum Fractional Thickness-weighted Barotropic Anomaly Zonal Acceleration', &

1655 'm s-2', conversion=us%L_T2_to_m_s2)

1656 if (cs%id_hf_u_BT_accel_2d > 0) call safe_alloc_ptr(cs%ADp%diag_hfrac_u,isdb,iedb,jsd,jed,nz)

1657

1658 cs%id_hf_v_BT_accel_2d = register_diag_field('ocean_model', 'hf_v_BT_accel_2d', diag%axesCv1, time, &

1659 'Depth-sum Fractional Thickness-weighted Barotropic Anomaly Meridional Acceleration', &

1660 'm s-2', conversion=us%L_T2_to_m_s2)

1661 if (cs%id_hf_v_BT_accel_2d > 0) call safe_alloc_ptr(cs%ADp%diag_hfrac_v,isd,ied,jsdb,jedb,nz)

1662

1663 cs%id_h_u_BT_accel = register_diag_field('ocean_model', 'h_u_BT_accel', diag%axesCuL, time, &

1664 'Thickness Multiplied Barotropic Anomaly Zonal Acceleration', &

1665 'm2 s-2', conversion=gv%H_to_m*us%L_T2_to_m_s2)

1666 if (cs%id_h_u_BT_accel > 0) call safe_alloc_ptr(cs%ADp%diag_hu,isdb,iedb,jsd,jed,nz)

1667

1668 cs%id_h_v_BT_accel = register_diag_field('ocean_model', 'h_v_BT_accel', diag%axesCvL, time, &

1669 'Thickness Multiplied Barotropic Anomaly Meridional Acceleration', &

1670 'm2 s-2', conversion=gv%H_to_m*us%L_T2_to_m_s2)

1671 if (cs%id_h_v_BT_accel > 0) call safe_alloc_ptr(cs%ADp%diag_hv,isd,ied,jsdb,jedb,nz)

1672

1673 cs%id_intz_u_BT_accel_2d = register_diag_field('ocean_model', 'intz_u_BT_accel_2d', diag%axesCu1, time, &

1674 'Depth-integral of Barotropic Anomaly Zonal Acceleration', &

1675 'm2 s-2', conversion=gv%H_to_m*us%L_T2_to_m_s2)

1676 if (cs%id_intz_u_BT_accel_2d > 0) call safe_alloc_ptr(cs%ADp%diag_hu,isdb,iedb,jsd,jed,nz)

1677

1678 cs%id_intz_v_BT_accel_2d = register_diag_field('ocean_model', 'intz_v_BT_accel_2d', diag%axesCv1, time, &

1679 'Depth-integral of Barotropic Anomaly Meridional Acceleration', &

1680 'm2 s-2', conversion=gv%H_to_m*us%L_T2_to_m_s2)

1681 if (cs%id_intz_v_BT_accel_2d > 0) call safe_alloc_ptr(cs%ADp%diag_hv,isd,ied,jsdb,jedb,nz)

1682

1683 cs%id_PFu_visc_rem = register_diag_field('ocean_model', 'PFu_visc_rem', diag%axesCuL, time, &

1684 'Zonal Pressure Force Acceleration multiplied by the viscous remnant', &

1685 'm s-2', conversion=us%L_T2_to_m_s2)

1686 if (cs%id_PFu_visc_rem > 0) call safe_alloc_ptr(cs%ADp%visc_rem_u,isdb,iedb,jsd,jed,nz)

1687 cs%id_PFv_visc_rem = register_diag_field('ocean_model', 'PFv_visc_rem', diag%axesCvL, time, &

1688 'Meridional Pressure Force Acceleration multiplied by the viscous remnant', &

1689 'm s-2', conversion=us%L_T2_to_m_s2)

1690 if (cs%id_PFv_visc_rem > 0) call safe_alloc_ptr(cs%ADp%visc_rem_v,isd,ied,jsdb,jedb,nz)

1691

1692 cs%id_CAu_visc_rem = register_diag_field('ocean_model', 'CAu_visc_rem', diag%axesCuL, time, &

1693 'Zonal Coriolis and Advective Acceleration multiplied by the viscous remnant', &

1694 'm s-2', conversion=us%L_T2_to_m_s2)

1695 if (cs%id_CAu_visc_rem > 0) call safe_alloc_ptr(cs%ADp%visc_rem_u,isdb,iedb,jsd,jed,nz)

1696 cs%id_CAv_visc_rem = register_diag_field('ocean_model', 'CAv_visc_rem', diag%axesCvL, time, &

1697 'Meridional Coriolis and Advective Acceleration multiplied by the viscous remnant', &

1698 'm s-2', conversion=us%L_T2_to_m_s2)

1699 if (cs%id_CAv_visc_rem > 0) call safe_alloc_ptr(cs%ADp%visc_rem_v,isd,ied,jsdb,jedb,nz)

1700

1701 cs%id_u_BT_accel_visc_rem = register_diag_field('ocean_model', 'u_BT_accel_visc_rem', diag%axesCuL, time, &

1702 'Barotropic Anomaly Zonal Acceleration multiplied by the viscous remnant', &

1703 'm s-2', conversion=us%L_T2_to_m_s2)

1704 if (cs%id_u_BT_accel_visc_rem > 0) call safe_alloc_ptr(cs%ADp%visc_rem_u,isdb,iedb,jsd,jed,nz)

1705 cs%id_v_BT_accel_visc_rem = register_diag_field('ocean_model', 'v_BT_accel_visc_rem', diag%axesCvL, time, &

1706 'Barotropic Anomaly Meridional Acceleration multiplied by the viscous remnant', &

1707 'm s-2', conversion=us%L_T2_to_m_s2)

1708 if (cs%id_v_BT_accel_visc_rem > 0) call safe_alloc_ptr(cs%ADp%visc_rem_v,isd,ied,jsdb,jedb,nz)

1709

1710 id_clock_cor = cpu_clock_id('(Ocean Coriolis & mom advection)', grain=clock_module)

1711 id_clock_continuity = cpu_clock_id('(Ocean continuity equation)', grain=clock_module)

1712 id_clock_pres = cpu_clock_id('(Ocean pressure force)', grain=clock_module)

1713 id_clock_vertvisc = cpu_clock_id('(Ocean vertical viscosity)', grain=clock_module)

1714 id_clock_horvisc = cpu_clock_id('(Ocean horizontal viscosity)', grain=clock_module)

1715 id_clock_mom_update = cpu_clock_id('(Ocean momentum increments)', grain=clock_module)

1716 id_clock_pass = cpu_clock_id('(Ocean message passing)', grain=clock_module)

1717 id_clock_btcalc = cpu_clock_id('(Ocean barotropic mode calc)', grain=clock_module)

1718 id_clock_btstep = cpu_clock_id('(Ocean barotropic mode stepping)', grain=clock_module)

1719 id_clock_btforce = cpu_clock_id('(Ocean barotropic forcing calc)', grain=clock_module)

1720

1721end subroutine initialize_dyn_split_rk2b

1722

1723

1724!> Close the dyn_split_RK2b module

1725subroutine end_dyn_split_rk2b(CS)

1726 type(mom_dyn_split_rk2b_cs), pointer :: cs !< module control structure

1727

1728 call barotropic_end(cs%barotropic_CSp)

1729

1730 call vertvisc_end(cs%vertvisc_CSp)

1731 deallocate(cs%vertvisc_CSp)

1732

1733 call hor_visc_end(cs%hor_visc)

1734 if (cs%calculate_SAL) call sal_end(cs%SAL_CSp)

1735 if (cs%use_tides) call tidal_forcing_end(cs%tides_CSp)

1736 call coriolisadv_end(cs%CoriolisAdv)

1737

1738 dealloc_(cs%diffu) ; dealloc_(cs%diffv)

1739 dealloc_(cs%CAu) ; dealloc_(cs%CAv)

1740 dealloc_(cs%CAu_pred) ; dealloc_(cs%CAv_pred)

1741 dealloc_(cs%PFu) ; dealloc_(cs%PFv)

1742

1743 if (associated(cs%taux_bot)) deallocate(cs%taux_bot)

1744 if (associated(cs%tauy_bot)) deallocate(cs%tauy_bot)

1745 dealloc_(cs%uhbt) ; dealloc_(cs%vhbt)

1746 dealloc_(cs%du_av_inst) ; dealloc_(cs%dv_av_inst)

1747 dealloc_(cs%u_accel_bt) ; dealloc_(cs%v_accel_bt)

1748 dealloc_(cs%visc_rem_u) ; dealloc_(cs%visc_rem_v)

1749

1750 dealloc_(cs%eta) ; dealloc_(cs%eta_PF) ; dealloc_(cs%pbce)

1751

1752 call dealloc_bt_cont_type(cs%BT_cont)

1753 deallocate(cs%AD_pred)

1754

1755 deallocate(cs)

1756end subroutine end_dyn_split_rk2b

1757

1758

1759!> \namespace mom_dynamics_split_rk2b

1760!!

1761!! This file time steps the adiabatic dynamic core by splitting

1762!! between baroclinic and barotropic modes. It uses a pseudo-second order

1763!! Runge-Kutta time stepping scheme for the baroclinic momentum

1764!! equation and a forward-backward coupling between the baroclinic

1765!! momentum and continuity equations. This split time-stepping

1766!! scheme is described in detail in Hallberg (JCP, 1997). Additional

1767!! issues related to exact tracer conservation and how to

1768!! ensure consistency between the barotropic and layered estimates

1769!! of the free surface height are described in Hallberg and

1770!! Adcroft (Ocean Modelling, 2009). This was the time stepping code

1771!! that is used for most GOLD applications, including GFDL's ESM2G

1772!! Earth system model, and all of the examples provided with the

1773!! MOM code (although several of these solutions are routinely

1774!! verified by comparison with the slower unsplit schemes).

1775!!

1776!! The subroutine step_MOM_dyn_split_RK2b actually does the time

1777!! stepping, while register_restarts_dyn_split_RK2b sets the fields

1778!! that are found in a full restart file with this scheme, and

1779!! initialize_dyn_split_RK2b initializes the cpu clocks that are

1780!! used in this module. For largely historical reasons, this module

1781!! does not have its own control structure, but shares the same

1782!! control structure with MOM.F90 and the other MOM_dynamics_...

1783!! modules.

1784

1785end module mom_dynamics_split_rk2b