src/core/MOM_dynamics_split_RK2.F90

portedportable, not yet portedexecuted, not portableexecutable, not hit by this run

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.
6		module MOM_dynamics_split_RK2
7
8		use MOM_variables, only : vertvisc_type, thermo_var_ptrs, porous_barrier_type
9		use MOM_variables, only : BT_cont_type, alloc_bt_cont_type, dealloc_bt_cont_type
10		use MOM_variables, only : accel_diag_ptrs, ocean_internal_state, cont_diag_ptrs
11		use MOM_forcing_type, only : mech_forcing
12
13		use MOM_checksum_packages, only : MOM_thermo_chksum, MOM_state_chksum, MOM_accel_chksum
14		use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
15		use MOM_cpu_clock, only : CLOCK_COMPONENT, CLOCK_SUBCOMPONENT
16		use MOM_cpu_clock, only : CLOCK_MODULE_DRIVER, CLOCK_MODULE, CLOCK_ROUTINE
17		use MOM_diabatic_driver, only : diabatic_CS, extract_diabatic_member
18		use MOM_diag_mediator, only : diag_mediator_init, enable_averages
19		use MOM_diag_mediator, only : disable_averaging, post_data, safe_alloc_ptr
20		use MOM_diag_mediator, only : post_product_u, post_product_sum_u
21		use MOM_diag_mediator, only : post_product_v, post_product_sum_v
22		use MOM_diag_mediator, only : register_diag_field, register_static_field
23		use MOM_diag_mediator, only : set_diag_mediator_grid, diag_ctrl, diag_update_remap_grids
24		use MOM_domains, only : To_South, To_West, To_All, CGRID_NE, SCALAR_PAIR
25		use MOM_domains, only : To_North, To_East, Omit_Corners
26		use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type
27		use MOM_domains, only : start_group_pass, complete_group_pass, pass_var, pass_vector
28		use MOM_debugging, only : hchksum, uvchksum, query_debugging_checks
29		use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, WARNING, is_root_pe
30		use MOM_error_handler, only : MOM_set_verbosity, callTree_showQuery
31		use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint
32		use MOM_file_parser, only : get_param, log_version, param_file_type
33		use MOM_get_input, only : directories
34		use MOM_io, only : vardesc, var_desc, EAST_FACE, NORTH_FACE
35		use MOM_restart, only : register_restart_field, register_restart_pair
36		use MOM_restart, only : query_initialized, set_initialized, save_restart
37		use MOM_restart, only : only_read_from_restarts
38		use MOM_restart, only : restart_init, is_new_run, MOM_restart_CS
39		use MOM_time_manager, only : time_type, operator(+)
40		use MOM_time_manager, only : operator(-), operator(>), operator(*), operator(/)
41
42		use MOM_ALE, only : ALE_CS, ALE_remap_velocities
43		use MOM_barotropic, only : barotropic_init, btstep, btcalc, bt_mass_source
44		use MOM_barotropic, only : register_barotropic_restarts, set_dtbt, barotropic_CS
45		use MOM_barotropic, only : barotropic_end
46		use MOM_boundary_update, only : update_OBC_data, update_OBC_CS
47		use MOM_continuity, only : continuity, continuity_CS
48		use MOM_continuity, only : continuity_init, continuity_stencil
49		use MOM_CoriolisAdv, only : CorAdCalc, CoriolisAdv_CS
50		use MOM_CoriolisAdv, only : CoriolisAdv_init, CoriolisAdv_end, CoriolisAdv_stencil
51		use MOM_CVMix_KPP, only : KPP_get_BLD, KPP_CS
52		use MOM_debugging, only : check_redundant
53		use MOM_energetic_PBL, only : energetic_PBL_get_MLD, energetic_PBL_CS
54		use MOM_grid, only : ocean_grid_type
55		use MOM_harmonic_analysis, only : HA_init, harmonic_analysis_CS
56		use MOM_hor_index, only : hor_index_type
57		use MOM_hor_visc, only : horizontal_viscosity, hor_visc_nkblock, hor_visc_CS, hor_visc_vel_stencil
58		use MOM_hor_visc, only : hor_visc_init, hor_visc_end
59		use MOM_interface_heights, only : thickness_to_dz, find_col_avg_SpV
60		use MOM_lateral_mixing_coeffs, only : VarMix_CS
61		use MOM_MEKE_types, only : MEKE_type
62		use MOM_open_boundary, only : ocean_OBC_type, radiation_open_bdry_conds
63		use MOM_open_boundary, only : open_boundary_zero_normal_flow, open_boundary_query
64		use MOM_open_boundary, only : open_boundary_test_extern_h, update_OBC_ramp
65		use MOM_open_boundary, only : copy_thickness_reservoirs
66		use MOM_open_boundary, only : update_segment_thickness_reservoirs
67		use MOM_PressureForce, only : PressureForce, PressureForce_CS
68		use MOM_PressureForce, only : PressureForce_init
69		use MOM_set_visc, only : set_viscous_ML, set_visc_CS
70		use MOM_stochastics, only : stochastic_CS
71		use MOM_thickness_diffuse, only : thickness_diffuse_CS
72		use MOM_self_attr_load, only : SAL_CS
73		use MOM_self_attr_load, only : SAL_init, SAL_end
74		use MOM_tidal_forcing, only : tidal_forcing_CS
75		use MOM_tidal_forcing, only : tidal_forcing_init, tidal_forcing_end
76		use MOM_unit_scaling, only : unit_scale_type
77		use MOM_vert_friction, only : vertvisc, vertvisc_coef, vertvisc_remnant
78		use MOM_vert_friction, only : vertvisc_init, vertvisc_end, vertvisc_CS
79		use MOM_vert_friction, only : updateCFLtruncationValue, vertFPmix
80		use MOM_verticalGrid, only : verticalGrid_type, get_thickness_units
81		use MOM_verticalGrid, only : get_flux_units, get_tr_flux_units
82		use MOM_wave_interface, only : wave_parameters_CS, Stokes_PGF
83
84		implicit none ; private
85
86		#include <MOM_memory.h>
87
88		!> MOM_dynamics_split_RK2 module control structure
89		type, public :: MOM_dyn_split_RK2_CS ; private
90		real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: &
91		CAu, & !< CAu = f*v - u.grad(u) [L T-2 ~> m s-2]
92		CAu_pred, & !< The predictor step value of CAu = f*v - u.grad(u) [L T-2 ~> m s-2]
93		PFu, & !< PFu = -dM/dx [L T-2 ~> m s-2]
94		PFu_Stokes, & !< PFu_Stokes = -d/dx int_r (u_L*duS/dr) [L T-2 ~> m s-2]
95		diffu !< Zonal acceleration due to convergence of the along-isopycnal stress tensor [L T-2 ~> m s-2]
96
97		real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: &
98		CAv, & !< CAv = -f*u - u.grad(v) [L T-2 ~> m s-2]
99		CAv_pred, & !< The predictor step value of CAv = -f*u - u.grad(v) [L T-2 ~> m s-2]
100		PFv, & !< PFv = -dM/dy [L T-2 ~> m s-2]
101		PFv_Stokes, & !< PFv_Stokes = -d/dy int_r (v_L*dvS/dr) [L T-2 ~> m s-2]
102		diffv !< Meridional acceleration due to convergence of the along-isopycnal stress tensor [L T-2 ~> m s-2]
103
104		real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: visc_rem_u
105		!< Both the fraction of the zonal momentum originally in a
106		!! layer that remains after a time-step of viscosity, and the
107		!! fraction of a time-step worth of a barotropic acceleration
108		!! that a layer experiences after viscosity is applied [nondim].
109		!! Nondimensional between 0 (at the bottom) and 1 (far above).
110		real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: u_accel_bt
111		!< The zonal layer accelerations due to the difference between
112		!! the barotropic accelerations and the baroclinic accelerations
113		!! that were fed into the barotopic calculation [L T-2 ~> m s-2]
114		real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: visc_rem_v
115		!< Both the fraction of the meridional momentum originally in
116		!! a layer that remains after a time-step of viscosity, and the
117		!! fraction of a time-step worth of a barotropic acceleration
118		!! that a layer experiences after viscosity is applied [nondim].
119		!! Nondimensional between 0 (at the bottom) and 1 (far above).
120		real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: v_accel_bt
121		!< The meridional layer accelerations due to the difference between
122		!! the barotropic accelerations and the baroclinic accelerations
123		!! that were fed into the barotopic calculation [L T-2 ~> m s-2]
124
125		! The following variables are only used with the split time stepping scheme.
126		real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: eta !< Instantaneous free surface height (in Boussinesq
127		!! mode) or column mass anomaly (in non-Boussinesq
128		!! mode) [H ~> m or kg m-2]
129		real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: u_av !< layer x-velocity with vertical mean replaced by
130		!! time-mean barotropic velocity over a baroclinic
131		!! timestep [L T-1 ~> m s-1]
132		real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: v_av !< layer y-velocity with vertical mean replaced by
133		!! time-mean barotropic velocity over a baroclinic
134		!! timestep [L T-1 ~> m s-1]
135		real ALLOCABLE_, dimension(NIMEM_,NJMEM_,NKMEM_) :: h_av !< arithmetic mean of two successive layer
136		!! thicknesses [H ~> m or kg m-2]
137		real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: eta_PF !< instantaneous SSH used in calculating PFu and
138		!! PFv [H ~> m or kg m-2]
139		real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_) :: uhbt !< average x-volume or mass flux determined by the
140		!! barotropic solver [H L2 T-1 ~> m3 s-1 or kg s-1].
141		!! uhbt is roughly equal to the vertical sum of uh.
142		real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_) :: vhbt !< average y-volume or mass flux determined by the
143		!! barotropic solver [H L2 T-1 ~> m3 s-1 or kg s-1].
144		!! vhbt is roughly equal to vertical sum of vh.
145		real ALLOCABLE_, dimension(NIMEM_,NJMEM_,NKMEM_) :: pbce !< pbce times eta gives the baroclinic pressure
146		!! anomaly in each layer due to free surface height
147		!! anomalies [L2 H-1 T-2 ~> m s-2 or m4 kg-1 s-2].
148		type(KPP_CS), pointer :: KPP_CSp => NULL() !< KPP control structure needed to ge
149		type(energetic_PBL_CS), pointer :: energetic_PBL_CSp => NULL() !< ePBL control structure
150
151		real, pointer, dimension(:,:) :: taux_bot => NULL() !< frictional x-bottom stress from the ocean
152		!! to the seafloor [R L Z T-2 ~> Pa]
153		real, pointer, dimension(:,:) :: tauy_bot => NULL() !< frictional y-bottom stress from the ocean
154		!! to the seafloor [R L Z T-2 ~> Pa]
155		type(BT_cont_type), pointer :: BT_cont => NULL() !< A structure with elements that describe the
156		!! effective summed open face areas as a function
157		!! of barotropic flow.
158
159		! This is to allow the previous, velocity-based coupling with between the
160		! baroclinic and barotropic modes.
161		logical :: BT_use_layer_fluxes !< If true, use the summed layered fluxes plus
162		!! an adjustment due to a changed barotropic
163		!! velocity in the barotropic continuity equation.
164		logical :: BT_adj_corr_mass_src !< If true, recalculates the barotropic mass source after
165		!! predictor step. This should make little difference in the
166		!! deep ocean but appears to help for vanished layers.
167		logical :: split_bottom_stress !< If true, provide the bottom stress
168		!! calculated by the vertical viscosity to the
169		!! barotropic solver.
170		logical :: dtbt_use_bt_cont !< If true, use BT_cont to calculate DTBT.
171		logical :: store_CAu !< If true, store the Coriolis and advective accelerations at the
172		!! end of the timestep for use in the next predictor step.
173		logical :: CAu_pred_stored !< If true, the Coriolis and advective accelerations at the
174		!! end of the timestep have been stored for use in the next
175		!! predictor step. This is used to accomodate various generations
176		!! of restart files.
177		logical :: calculate_SAL !< If true, calculate self-attraction and loading.
178		logical :: use_tides !< If true, tidal forcing is enabled.
179		logical :: use_HA !< If true, perform inline harmonic analysis.
180		logical :: remap_aux !< If true, apply ALE remapping to all of the auxiliary 3-D
181		!! variables that are needed to reproduce across restarts,
182		!! similarly to what is done with the primary state variables.
183
184		real :: be !< A nondimensional number from 0.5 to 1 that controls
185		!! the backward weighting of the time stepping scheme [nondim]
186		real :: begw !< A nondimensional number from 0 to 1 that controls
187		!! the extent to which the treatment of gravity waves
188		!! is forward-backward (0) or simulated backward
189		!! Euler (1) [nondim]. 0 is often used.
190		real :: Cemp_NL !< Empirical coefficient of non-local momentum mixing [nondim]
191		logical :: debug !< If true, write verbose checksums for debugging purposes.
192		logical :: debug_OBC !< If true, do additional calls resetting values to help debug the correctness
193		!! of the open boundary condition code.
194		logical :: fpmix !< If true, add non-local momentum flux increments and diffuse down the Eulerian gradient.
195		logical :: module_is_initialized = .false. !< Record whether this module has been initialized.
196		logical :: visc_rem_dt_bug = .true. !< If true, recover a bug that uses dt_pred rather than dt for vertvisc_rem
197		!! at the end of predictor.
198
199		!>@{ Diagnostic IDs
200		integer :: id_uh = -1, id_vh = -1
201		integer :: id_umo = -1, id_vmo = -1
202		integer :: id_umo_2d = -1, id_vmo_2d = -1
203		integer :: id_PFu = -1, id_PFv = -1
204		integer :: id_CAu = -1, id_CAv = -1
205		integer :: id_ueffA = -1, id_veffA = -1
206		! integer :: id_hf_PFu = -1, id_hf_PFv = -1
207		integer :: id_h_PFu = -1, id_h_PFv = -1
208		integer :: id_hf_PFu_2d = -1, id_hf_PFv_2d = -1
209		integer :: id_intz_PFu_2d = -1, id_intz_PFv_2d = -1
210		integer :: id_PFu_visc_rem = -1, id_PFv_visc_rem = -1
211		! integer :: id_hf_CAu = -1, id_hf_CAv = -1
212		integer :: id_h_CAu = -1, id_h_CAv = -1
213		integer :: id_hf_CAu_2d = -1, id_hf_CAv_2d = -1
214		integer :: id_intz_CAu_2d = -1, id_intz_CAv_2d = -1
215		integer :: id_CAu_visc_rem = -1, id_CAv_visc_rem = -1
216		integer :: id_deta_dt = -1
217
218		! Split scheme only.
219		integer :: id_uav = -1, id_vav = -1
220		integer :: id_u_BT_accel = -1, id_v_BT_accel = -1
221		! integer :: id_hf_u_BT_accel = -1, id_hf_v_BT_accel = -1
222		integer :: id_h_u_BT_accel = -1, id_h_v_BT_accel = -1
223		integer :: id_hf_u_BT_accel_2d = -1, id_hf_v_BT_accel_2d = -1
224		integer :: id_intz_u_BT_accel_2d = -1, id_intz_v_BT_accel_2d = -1
225		integer :: id_u_BT_accel_visc_rem = -1, id_v_BT_accel_visc_rem = -1
226		!>@}
227
228		type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to regulate the
229		!! timing of diagnostic output.
230		type(accel_diag_ptrs), pointer :: ADp => NULL() !< A structure pointing to the various
231		!! accelerations in the momentum equations,
232		!! which can later be used to calculate
233		!! derived diagnostics like energy budgets.
234		type(accel_diag_ptrs), pointer :: AD_pred => NULL() !< A structure pointing to the various
235		!! predictor step accelerations in the momentum equations,
236		!! which can be used to debug truncations.
237		type(cont_diag_ptrs), pointer :: CDp => NULL() !< A structure with pointers to various
238		!! terms in the continuity equations,
239		!! which can later be used to calculate
240		!! derived diagnostics like energy budgets.
241
242		! The remainder of the structure points to child subroutines' control structures.
243		!> A pointer to the horizontal viscosity control structure
244		type(hor_visc_CS) :: hor_visc
245		!> A pointer to the continuity control structure
246		type(continuity_CS) :: continuity_CSp
247		!> The CoriolisAdv control structure
248		type(CoriolisAdv_CS) :: CoriolisAdv
249		!> A pointer to the PressureForce control structure
250		type(PressureForce_CS) :: PressureForce_CSp
251		!> A pointer to a structure containing interface height diffusivities
252		type(vertvisc_CS), pointer :: vertvisc_CSp => NULL()
253		!> A pointer to the set_visc control structure
254		type(set_visc_CS), pointer :: set_visc_CSp => NULL()
255		!> A pointer to the barotropic stepping control structure
256		type(barotropic_CS) :: barotropic_CSp
257		!> A pointer to the SAL control structure
258		type(SAL_CS) :: SAL_CSp
259		!> A pointer to the tidal forcing control structure
260		type(tidal_forcing_CS) :: tides_CSp
261		!> A pointer to the harmonic analysis control structure
262		type(harmonic_analysis_CS) :: HA_CSp
263		!> A pointer to the ALE control structure.
264		type(ALE_CS), pointer :: ALE_CSp => NULL()
265
266		type(ocean_OBC_type), pointer :: OBC => NULL() !< A pointer to an open boundary
267		!! condition type that specifies whether, where, and what open boundary
268		!! conditions are used. If no open BCs are used, this pointer stays
269		!! nullified. Flather OBCs use open boundary_CS as well.
270		!> A pointer to the update_OBC control structure
271		type(update_OBC_CS), pointer :: update_OBC_CSp => NULL()
272
273		type(group_pass_type) :: pass_eta !< Structure for group halo pass
274		type(group_pass_type) :: pass_visc_rem !< Structure for group halo pass
275		type(group_pass_type) :: pass_uvp !< Structure for group halo pass
276		type(group_pass_type) :: pass_hp_uv !< Structure for group halo pass
277		type(group_pass_type) :: pass_uv !< Structure for group halo pass
278		type(group_pass_type) :: pass_h !< Structure for group halo pass
279		type(group_pass_type) :: pass_av_uvh !< Structure for group halo pass
280
281		end type MOM_dyn_split_RK2_CS
282
283
284		public step_MOM_dyn_split_RK2
285		public register_restarts_dyn_split_RK2
286		public initialize_dyn_split_RK2
287		public remap_dyn_split_RK2_aux_vars
288		public init_dyn_split_RK2_diabatic
289		public end_dyn_split_RK2
290
291		!>@{ CPU time clock IDs
292		integer :: id_clock_Cor, id_clock_pres, id_clock_vertvisc
293		integer :: id_clock_horvisc, id_clock_mom_update
294		integer :: id_clock_continuity, id_clock_thick_diff
295		integer :: id_clock_btstep, id_clock_btcalc, id_clock_btforce
296		integer :: id_clock_pass, id_clock_pass_init
297		!>@}
298
299		contains
300
301		!> RK2 splitting for time stepping MOM adiabatic dynamics
302	24	subroutine step_MOM_dyn_split_RK2(u_inst, v_inst, h, tv, visc, Time_local, dt, forces, &
303	24	p_surf_begin, p_surf_end, uh, vh, uhtr, vhtr, eta_av, G, GV, US, CS, &
304		calc_dtbt, VarMix, MEKE, thickness_diffuse_CSp, pbv, STOCH, Waves)
305		type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure
306		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
307		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
308		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
309		target, intent(inout) :: u_inst !< Instantaneous zonal velocity [L T-1 ~> m s-1]
310		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
311		target, intent(inout) :: v_inst !< Instantaneous meridional velocity [L T-1 ~> m s-1]
312		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
313		intent(inout) :: h !< Layer thickness [H ~> m or kg m-2]
314		type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic type
315		type(vertvisc_type), intent(inout) :: visc !< Vertical visc, bottom drag, and related
316		type(time_type), intent(in) :: Time_local !< Model time at end of time step
317		real, intent(in) :: dt !< Baroclinic dynamics time step [T ~> s]
318		type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces
319		real, dimension(:,:), pointer :: p_surf_begin !< Surface pressure at the start of this dynamic
320		!! time step [R L2 T-2 ~> Pa]
321		real, dimension(:,:), pointer :: p_surf_end !< Surface pressure at the end of this dynamic
322		!! time step [R L2 T-2 ~> Pa]
323		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
324		target, intent(inout) :: uh !< Zonal volume or mass transport
325		!! [H L2 T-1 ~> m3 s-1 or kg s-1]
326		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
327		target, intent(inout) :: vh !< Meridional volume or mass transport
328		!! [H L2 T-1 ~> m3 s-1 or kg s-1]
329		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
330		intent(inout) :: uhtr !< Accumulated zonal volume or mass transport
331		!! since last tracer advection [H L2 ~> m3 or kg]
332		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
333		intent(inout) :: vhtr !< Accumulated meridional volume or mass transport
334		!! since last tracer advection [H L2 ~> m3 or kg]
335		real, dimension(SZI_(G),SZJ_(G)), intent(out) :: eta_av !< Free surface height or column mass
336		!! averaged over time step [H ~> m or kg m-2]
337		type(MOM_dyn_split_RK2_CS), pointer :: CS !< Module control structure
338		logical, intent(in) :: calc_dtbt !< If true, recalculate the barotropic time step
339		type(VarMix_CS), intent(inout) :: VarMix !< Variable mixing control structure
340		type(MEKE_type), intent(inout) :: MEKE !< MEKE fields
341		type(thickness_diffuse_CS), intent(inout) :: thickness_diffuse_CSp !< Pointer to a structure containing
342		!! interface height diffusivities
343		type(porous_barrier_type), intent(in) :: pbv !< porous barrier fractional cell metrics
344		type(stochastic_CS), optional, intent(inout) :: STOCH !< Stochastic control structure
345		type(wave_parameters_CS), optional, pointer :: Waves !< A pointer to a structure containing
346		!! fields related to the surface wave conditions
347
348		! local variables
349	48	real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: up ! Predicted zonal velocity [L T-1 ~> m s-1].
350	48	real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: vp ! Predicted meridional velocity [L T-1 ~> m s-1].
351	48	real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: hp ! Predicted thickness [H ~> m or kg m-2].
352	48	real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: dz ! Distance between the interfaces around a layer [Z ~> m]
353	48	real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: ueffA ! Effective Area of U-Faces [H L ~> m2]
354	48	real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: veffA ! Effective Area of V-Faces [H L ~> m2]
355	48	real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: u_bc_accel ! The summed zonal baroclinic accelerations
356		! of each layer calculated by the non-barotropic
357		! part of the model [L T-2 ~> m s-2]
358	48	real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: v_bc_accel ! The summed meridional baroclinic accelerations
359		! of each layer calculated by the non-barotropic
360		! part of the model [L T-2 ~> m s-2]
361
362	48	real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), target :: uh_in ! The zonal mass transports that would be
363		! obtained using the initial velocities [H L2 T-1 ~> m3 s-1 or kg s-1]
364	48	real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), target :: vh_in ! The meridional mass transports that would be
365		! obtained using the initial velocities [H L2 T-1 ~> m3 s-1 or kg s-1]
366
367	48	real, dimension(SZI_(G),SZJ_(G)) :: eta_pred ! The predictor value of the free surface height
368		! or column mass [H ~> m or kg m-2]
369	48	real, dimension(SZI_(G),SZJ_(G)) :: SpV_avg ! The column averaged specific volume [R-1 ~> m3 kg-1]
370	48	real, dimension(SZI_(G),SZJ_(G)) :: deta_dt ! A diagnostic of the time derivative of the free surface
371		! height or column mass [H T-1 ~> m s-1 or kg m-2 s-1]
372
373	48	real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: u_old_rad_OBC ! The starting zonal velocities, which are
374		! saved for use in the radiation open boundary condition code [L T-1 ~> m s-1]
375	48	real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: v_old_rad_OBC ! The starting meridional velocities, which are
376		! saved for use in the radiation open boundary condition code [L T-1 ~> m s-1]
377
378		! GMM, TODO: make these allocatable?
379	48	real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: uold ! u-velocity before vert_visc is applied, for fpmix
380		! [L T-1 ~> m s-1]
381	48	real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: vold ! v-velocity before vert_visc is applied, for fpmix
382		! [L T-1 ~> m s-1]
383		real :: pres_to_eta ! A factor that converts pressures to the units of eta
384		! [H T2 R-1 L-2 ~> m Pa-1 or kg m-2 Pa-1]
385		real, pointer, dimension(:,:) :: &
386		p_surf => NULL(), & ! A pointer to the surface pressure [R L2 T-2 ~> Pa]
387		eta_PF_start => NULL(), & ! The value of eta that corresponds to the starting pressure
388		! for the barotropic solver [H ~> m or kg m-2]
389		taux_bot => NULL(), & ! A pointer to the zonal bottom stress in some cases [R L Z T-2 ~> Pa]
390		tauy_bot => NULL(), & ! A pointer to the meridional bottom stress in some cases [R L Z T-2 ~> Pa]
391		! This pointer is just used as shorthand for CS%eta.
392		eta => NULL() ! A pointer to the instantaneous free surface height (in Boussinesq
393		! mode) or column mass anomaly (in non-Boussinesq mode) [H ~> m or kg m-2]
394
395		real, pointer, dimension(:,:,:) :: &
396		! These pointers are used to alter which fields are passed to btstep with various options:
397		u_ptr => NULL(), & ! A pointer to a zonal velocity [L T-1 ~> m s-1]
398		v_ptr => NULL(), & ! A pointer to a meridional velocity [L T-1 ~> m s-1]
399		uh_ptr => NULL(), & ! A pointer to a zonal volume or mass transport [H L2 T-1 ~> m3 s-1 or kg s-1]
400		vh_ptr => NULL(), & ! A pointer to a meridional volume or mass transport [H L2 T-1 ~> m3 s-1 or kg s-1]
401		! These pointers are just used as shorthand for CS%u_av, CS%v_av, and CS%h_av.
402	24	u_av, & ! The zonal velocity time-averaged over a time step [L T-1 ~> m s-1].
403	24	v_av, & ! The meridional velocity time-averaged over a time step [L T-1 ~> m s-1].
404	24	h_av ! The layer thickness time-averaged over a time step [H ~> m or kg m-2].
405
406	48	real, dimension(SZI_(G),SZJ_(G)) :: hbl ! Boundary layer depth from Cvmix [H ~> m or kg m-2]
407		real :: dt_pred ! The time step for the predictor part of the baroclinic time stepping [T ~> s].
408		real :: Idt_bc ! Inverse of the baroclinic timestep [T-1 ~> s-1]
409		logical :: dyn_p_surf
410		logical :: debug_redundant ! If true, check redundant values on PE boundaries when debugging
411		logical :: BT_cont_BT_thick ! If true, use the BT_cont_type to estimate the
412		! relative weightings of the layers in calculating
413		! the barotropic accelerations.
414		logical :: Use_Stokes_PGF ! If true, add Stokes PGF to hydrostatic PGF
415		!---For group halo pass
416		logical :: showCallTree, sym
417		logical :: lFPpost ! Used to only post diagnostics in vertFPmix when fpmix=true and
418		! in the corrector step (not the predict)
419		integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
420		integer :: cont_stencil, obc_stencil, vel_stencil
421	48	real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_tmp ! temporary copy of Layer thickness [H ~> m or kg m-2]
422		integer :: cor_stencil
423
424	24	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
425	24	Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
426	24	u_av => CS%u_av ; v_av => CS%v_av ; h_av => CS%h_av ; eta => CS%eta
427
428	24	Idt_bc = 1.0 / dt
429
430	24	sym = G%Domain%symmetric ! switch to include symmetric domain in checksums
431
432	24	showCallTree = callTree_showQuery()
433	0	if (showCallTree) call callTree_enter("step_MOM_dyn_split_RK2(), MOM_dynamics_split_RK2.F90")
434
435		!$omp target enter data map(alloc: u_bc_accel, v_bc_accel, eta_pred, uh_in, vh_in)
436		!$omp target enter data map(alloc: up, vp, hp, dz, h_tmp)
437
438	3120	do concurrent (k=1:nz, j=G%jsd:G%jed, I=G%IsdB:G%IedB)
439	16003248	up(I,j,k) = 0.0
440		enddo
441	3096	do concurrent (k=1:nz, J=G%JsdB:G%JedB, i=G%isd:G%ied)
442	16112664	vp(i,J,k) = 0.0
443		enddo
444	3096	do concurrent (k=1:nz, j=G%jsd:G%jed, i=G%isd:G%ied)
445	15879192	hp(i,j,k) = h(i,j,k)
446		enddo
447		! TODO: hp needs accurate +/-2 halos.
448
449		! Update CFL truncation value as function of time
450	24	call updateCFLtruncationValue(Time_local, CS%vertvisc_CSp, US)
451
452	24	if (CS%debug) then
453	0	call query_debugging_checks(do_redundant=debug_redundant)
454	0	call MOM_state_chksum("Start predictor ", u_inst, v_inst, h, uh, vh, G, GV, US, symmetric=sym)
455	0	if (debug_redundant) then
456	0	call check_redundant("Start predictor u ", u_inst, v_inst, G, unscale=US%L_T_to_m_s)
457	0	call check_redundant("Start predictor uh ", uh, vh, G, unscale=GV%H_to_MKSUS%L_to_m2US%s_to_T)
458		endif
459		endif
460
461	24	dyn_p_surf = associated(p_surf_begin) .and. associated(p_surf_end)
462	24	if (dyn_p_surf) then
463	0	p_surf => p_surf_end
464	0	call safe_alloc_ptr(eta_PF_start,G%isd,G%ied,G%jsd,G%jed)
465	0	eta_PF_start(:,:) = 0.0
466		else
467	24	p_surf => forces%p_surf
468		endif
469		! TODO: This should probably be resolved in step_MOM, if not higher. But
470		! p_surf setup is a bit complicated, and may even depend on the driver, so
471		! for now we alloc/delete in the dynamic core step.
472		!$omp target enter data map(to: p_surf) if (associated(p_surf))
473
474	24	if (associated(CS%OBC)) then
475	0	if (CS%debug_OBC) call open_boundary_test_extern_h(G, GV, CS%OBC, h)
476
477		! Update OBC ramp value as function of time
478	0	call update_OBC_ramp(Time_local, CS%OBC, US)
479
480	0	do k=1,nz ; do j=G%jsd,G%jed ; do I=G%IsdB,G%IedB
481	0	u_old_rad_OBC(I,j,k) = u_av(I,j,k)
482		enddo ; enddo ; enddo
483	0	do k=1,nz ; do J=G%JsdB,G%JedB ; do i=G%isd,G%ied
484	0	v_old_rad_OBC(i,J,k) = v_av(i,J,k)
485		enddo ; enddo ; enddo
486		endif
487
488	24	BT_cont_BT_thick = .false.
489	24	if (associated(CS%BT_cont)) BT_cont_BT_thick = &
490	24	(allocated(CS%BT_cont%h_u) .and. allocated(CS%BT_cont%h_v))
491
492	24	if (CS%split_bottom_stress) then
493	0	taux_bot => CS%taux_bot ; tauy_bot => CS%tauy_bot
494		endif
495
496		!--- begin set up for group halo pass
497
498	24	cont_stencil = continuity_stencil(CS%continuity_CSp)
499	24	obc_stencil = 2
500	24	cor_stencil = CoriolisAdv_stencil(CS%CoriolisAdv)
501	24	if (associated(CS%OBC)) then
502	0	if (CS%OBC%oblique_BCs_exist_globally) obc_stencil = 3
503		endif
504	24	vel_stencil = max(2, obc_stencil, hor_visc_vel_stencil(CS%hor_visc))
505	24	call cpu_clock_begin(id_clock_pass)
506	24	call create_group_pass(CS%pass_eta, eta, G%Domain, halo=1)
507		call create_group_pass(CS%pass_visc_rem, CS%visc_rem_u, CS%visc_rem_v, G%Domain, &
508	24	To_All+SCALAR_PAIR, CGRID_NE, halo=max(1,cont_stencil))
509	24	call create_group_pass(CS%pass_uvp, up, vp, G%Domain, halo=max(1,cont_stencil))
510	24	call create_group_pass(CS%pass_uv, u_inst, v_inst, G%Domain, halo=max(2,cont_stencil))
511
512	24	call create_group_pass(CS%pass_hp_uv, hp, G%Domain, halo=cor_stencil)
513	24	call create_group_pass(CS%pass_hp_uv, u_av, v_av, G%Domain, halo=max(cor_stencil,vel_stencil))
514	24	call create_group_pass(CS%pass_hp_uv, uh(:,:,:), vh(:,:,:), G%Domain, halo=max(cor_stencil,vel_stencil))
515	24	call create_group_pass(CS%pass_h, h, g%domain, halo=max(cor_stencil,cont_stencil))
516	24	call create_group_pass(CS%pass_av_uvh, u_av, v_av, g%domain, halo=max(cor_stencil,vel_stencil))
517	24	call create_group_pass(CS%pass_av_uvh, uh(:,:,:), vh(:,:,:), G%Domain, halo=max(cor_stencil,vel_stencil))
518
519	24	call cpu_clock_end(id_clock_pass)
520		!--- end set up for group halo pass
521
522		! PFu = d/dx M(h,T,S)
523		! pbce = dM/deta
524	24	if (CS%begw == 0.0) call enable_averages(dt, Time_local, CS%diag)
525	24	call cpu_clock_begin(id_clock_pres)
526
527		call PressureForce(h, tv, CS%PFu, CS%PFv, G, GV, US, CS%PressureForce_CSp, &
528	24	CS%ALE_CSp, CS%ADp, p_surf, CS%pbce, CS%eta_PF)
529
530	24	if (dyn_p_surf) then
531		!$omp target update from(CS%eta_PF)
532	0	pres_to_eta = 1.0 / (GV%g_Earth * GV%H_to_RZ)
533		!$OMP parallel do default(shared)
534	0	do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
535	0	eta_PF_start(i,j) = CS%eta_PF(i,j) - pres_to_eta * (p_surf_begin(i,j) - p_surf_end(i,j))
536		enddo ; enddo
537		endif
538		! Stokes shear force contribution to pressure gradient
539	24	Use_Stokes_PGF = present(Waves)
540	24	if (Use_Stokes_PGF) then
541	24	Use_Stokes_PGF = associated(Waves)
542	24	if (Use_Stokes_PGF) Use_Stokes_PGF = Waves%Stokes_PGF
543	24	if (Use_Stokes_PGF) then
544		!$omp target update from(h, u_inst, v_inst, CS%PFu, CS%PFv)
545	0	call thickness_to_dz(h, tv, dz, G, GV, US, halo_size=1)
546	0	call Stokes_PGF(G, GV, US, dz, u_inst, v_inst, CS%PFu_Stokes, CS%PFv_Stokes, Waves)
547
548		! We are adding Stokes_PGF to hydrostatic PGF here. The diag PFu/PFv
549		! will therefore report the sum total PGF and we avoid other
550		! modifications in the code. The PFu_Stokes is output within the waves routines.
551	0	if (.not.Waves%Passive_Stokes_PGF) then
552	0	do k=1,nz
553	0	do j=js,je ; do I=Isq,Ieq
554	0	CS%PFu(I,j,k) = CS%PFu(I,j,k) + CS%PFu_Stokes(I,j,k)
555		enddo ; enddo
556		enddo
557	0	do k=1,nz
558	0	do J=Jsq,Jeq ; do i=is,ie
559	0	CS%PFv(i,J,k) = CS%PFv(i,J,k) + CS%PFv_Stokes(i,J,k)
560		enddo ; enddo
561		enddo
562		endif
563		endif
564		endif
565	24	call cpu_clock_end(id_clock_pres)
566	24	call disable_averaging(CS%diag)
567	24	if (showCallTree) call callTree_wayPoint("done with PressureForce (step_MOM_dyn_split_RK2)")
568
569	24	if (associated(CS%OBC)) then ; if (CS%OBC%update_OBC) then
570	0	call update_OBC_data(CS%OBC, G, GV, US, tv, h, CS%update_OBC_CSp, Time_local)
571		endif ; endif
572
573	24	if (associated(CS%OBC) .and. CS%debug_OBC) then
574		!$omp target update from (CS%PFu, CS%PFv)
575	0	call open_boundary_zero_normal_flow(CS%OBC, G, GV, CS%PFu, CS%PFv)
576		endif
577
578	24	if (G%nonblocking_updates) then
579		!$omp target update from(eta)
580	0	call start_group_pass(CS%pass_eta, G%Domain, clock=id_clock_pass)
581		endif
582
583		! CAu = -(f+zeta_av)/h_av vh + d/dx KE_av
584	24	if (.not.CS%CAu_pred_stored) then
585		! Calculate a predictor-step estimate of the Coriolis and momentum advection terms,
586		! if it was not already stored from the end of the previous time step.
587	0	call cpu_clock_begin(id_clock_Cor)
588		!$omp target update to(u_av, v_av, h_av, uh, vh)
589		call CorAdCalc(u_av, v_av, h_av, uh, vh, CS%CAu_pred, CS%CAv_pred, CS%OBC, CS%AD_pred, &
590	0	G, GV, US, CS%CoriolisAdv, pbv, Waves=Waves)
591	0	call cpu_clock_end(id_clock_Cor)
592	0	if (showCallTree) call callTree_wayPoint("done with CorAdCalc (step_MOM_dyn_split_RK2)")
593		endif
594
595		! u_bc_accel = CAu + PFu + diffu(u[n-1])
596	24	call cpu_clock_begin(id_clock_btforce)
597
598	24	do concurrent (k=1:nz, j=js:je, I=Isq:Ieq)
599	13245168	u_bc_accel(I,j,k) = (CS%CAu_pred(I,j,k) + CS%PFu(I,j,k)) + CS%diffu(I,j,k)
600		enddo
601	24	do concurrent (k=1:nz, J=Jsq:Jeq, i=is:ie)
602	13354584	v_bc_accel(i,J,k) = (CS%CAv_pred(i,J,k) + CS%PFv(i,J,k)) + CS%diffv(i,J,k)
603		enddo
604
605	24	if (associated(CS%OBC)) then
606		!$omp target update from(u_bc_accel, v_bc_accel)
607	0	call open_boundary_zero_normal_flow(CS%OBC, G, GV, u_bc_accel, v_bc_accel)
608		!$omp target update to(u_bc_accel, v_bc_accel)
609		endif
610	24	call cpu_clock_end(id_clock_btforce)
611
612	24	if (CS%debug) then
613		!$omp target update from(CS%CAu_pred, CS%CAv_pred)
614		!$omp target update from(CS%PFu, CS%PFv, CS%pbce)
615		!$omp target update from(CS%diffu, CS%diffv)
616		!$omp target update from(u_bc_accel, v_bc_accel)
617		call MOM_accel_chksum("pre-btstep accel", CS%CAu_pred, CS%CAv_pred, CS%PFu, CS%PFv, &
618		CS%diffu, CS%diffv, G, GV, US, CS%pbce, u_bc_accel, v_bc_accel, &
619	0	symmetric=sym)
620	0	if (debug_redundant) then
621	0	call check_redundant("pre-btstep CS%CA ", CS%CAu_pred, CS%CAv_pred, G, unscale=US%L_T2_to_m_s2)
622	0	call check_redundant("pre-btstep CS%PF ", CS%PFu, CS%PFv, G, unscale=US%L_T2_to_m_s2)
623	0	call check_redundant("pre-btstep CS%diff ", CS%diffu, CS%diffv, G, unscale=US%L_T2_to_m_s2)
624	0	call check_redundant("pre-btstep u_bc_accel ", u_bc_accel, v_bc_accel, G, unscale=US%L_T2_to_m_s2)
625		endif
626		endif
627
628	24	call cpu_clock_begin(id_clock_vertvisc)
629
630	24	do concurrent (k=1:nz, j=js:je, I=Isq:Ieq)
631	13245168	up(I,j,k) = G%mask2dCu(I,j) * (u_inst(I,j,k) + dt * u_bc_accel(I,j,k))
632		enddo
633	24	do concurrent (k=1:nz, J=Jsq:Jeq, i=is:ie)
634	13354584	vp(i,J,k) = G%mask2dCv(i,J) * (v_inst(i,J,k) + dt * v_bc_accel(i,J,k))
635		enddo
636
637	24	call enable_averages(dt, Time_local, CS%diag)
638		! NOTE: this is on CPU and conditionally called (using `if (..) return`)
639		! It contains GPU/CPU data transfers for [uv]_inst and visc fields.
640	24	call set_viscous_ML(u_inst, v_inst, h, tv, forces, visc, dt, G, GV, US, CS%set_visc_CSp)
641	24	call disable_averaging(CS%diag)
642
643	24	if (CS%debug) then
644		!$omp target update from(up, vp)
645	0	call uvchksum("before vertvisc: up", up, vp, G%HI, haloshift=0, symmetric=sym, unscale=US%L_T_to_m_s)
646		endif
647
648	24	call thickness_to_dz(h, tv, dz, G, GV, US, halo_size=1, do_offload=.true.)
649
650	24	call vertvisc_coef(up, vp, h, dz, forces, visc, tv, dt, G, GV, US, CS%vertvisc_CSp, CS%OBC, VarMix)
651	24	call vertvisc_remnant(visc, CS%visc_rem_u, CS%visc_rem_v, dt, G, GV, US, CS%vertvisc_CSp)
652
653	24	call cpu_clock_end(id_clock_vertvisc)
654	24	if (showCallTree) call callTree_wayPoint("done with vertvisc_coef (step_MOM_dyn_split_RK2)")
655
656	24	call cpu_clock_begin(id_clock_pass)
657	24	if (G%nonblocking_updates) then
658	0	call complete_group_pass(CS%pass_eta, G%Domain)
659		!$omp target update to(eta)
660		!$omp target update from(CS%visc_rem_u, CS%visc_rem_v)
661	0	call start_group_pass(CS%pass_visc_rem, G%Domain)
662		else
663	24	call do_group_pass(CS%pass_eta, G%Domain, omp_offload=.true.)
664	24	call do_group_pass(CS%pass_visc_rem, G%Domain, omp_offload=.true.)
665		endif
666	24	call cpu_clock_end(id_clock_pass)
667
668	24	call cpu_clock_begin(id_clock_btcalc)
669		! Calculate the relative layer weights for determining barotropic quantities.
670	24	if (.not.BT_cont_BT_thick) then
671	0	call btcalc(h, G, GV, CS%barotropic_CSp, OBC=CS%OBC)
672		endif
673	24	call bt_mass_source(h, eta, .true., G, GV, CS%barotropic_CSp)
674
675	210552	SpV_avg(:,:) = 0.0
676	24	if ((.not.GV%Boussinesq) .and. associated(CS%OBC)) then
677		! Determine the column average specific volume if it is needed due to the
678		! use of Flather open boundary conditions in non-Boussinesq mode.
679	0	if (open_boundary_query(CS%OBC, apply_Flather_OBC=.true.)) &
680	0	call find_col_avg_SpV(h, SpV_avg, tv, G, GV, US)
681		endif
682	24	call cpu_clock_end(id_clock_btcalc)
683
684	24	if (G%nonblocking_updates) then
685	0	call complete_group_pass(CS%pass_visc_rem, G%Domain, clock=id_clock_pass)
686		!$omp target update to(CS%visc_rem_u, CS%visc_rem_v)
687		endif
688
689	24	if (associated(CS%OBC)) &
690	0	call copy_thickness_reservoirs(CS%OBC, G, GV)
691
692		! u_accel_bt = layer accelerations due to barotropic solver
693	24	if (associated(CS%BT_cont) .or. CS%BT_use_layer_fluxes) then
694	24	call cpu_clock_begin(id_clock_continuity)
695		call continuity(u_inst, v_inst, h, hp, uh_in, vh_in, dt, G, GV, US, CS%continuity_CSp, CS%OBC, pbv, &
696	24	visc_rem_u=CS%visc_rem_u, visc_rem_v=CS%visc_rem_v, BT_cont=CS%BT_cont)
697
698	24	call cpu_clock_end(id_clock_continuity)
699	24	if (BT_cont_BT_thick) then
700
701		call btcalc(h, G, GV, CS%barotropic_CSp, CS%BT_cont%h_u, CS%BT_cont%h_v, &
702	24	OBC=CS%OBC)
703
704		endif
705	24	if (showCallTree) call callTree_wayPoint("done with continuity[BT_cont] (step_MOM_dyn_split_RK2)")
706		endif
707
708	24	if (CS%BT_use_layer_fluxes) then
709	24	uh_ptr => uh_in ; vh_ptr => vh_in ; u_ptr => u_inst ; v_ptr => v_inst
710		endif
711
712	24	call cpu_clock_begin(id_clock_btstep)
713	24	if (calc_dtbt) then
714	13	if (CS%dtbt_use_bt_cont .and. associated(CS%BT_cont)) then
715	0	call set_dtbt(G, GV, US, CS%barotropic_CSp, CS%pbce, BT_cont=CS%BT_cont)
716		else
717		! In the following call, eta is only used when NONLINEAR_BT_CONTINUITY is True. Otherwise, dtbt is effectively
718		! calculated with eta=0. Note that NONLINEAR_BT_CONTINUITY is False if BT_CONT is used, which is the default.
719	13	call set_dtbt(G, GV, US, CS%barotropic_CSp, CS%pbce, eta=eta)
720		endif
721		endif
722
723	24	if (showCallTree) call callTree_enter("btstep(), MOM_barotropic.F90")
724		! This is the predictor step call to btstep.
725		! The CS%ADp argument here stores the weights for certain integrated diagnostics.
726		call btstep(u_inst, v_inst, eta, dt, u_bc_accel, v_bc_accel, forces, CS%pbce, CS%eta_PF, u_av, v_av, &
727		CS%u_accel_bt, CS%v_accel_bt, eta_pred, CS%uhbt, CS%vhbt, G, GV, US, &
728		CS%barotropic_CSp, CS%visc_rem_u, CS%visc_rem_v, SpV_avg, CS%ADp, CS%OBC, CS%BT_cont, &
729	24	eta_PF_start, taux_bot, tauy_bot, uh_ptr, vh_ptr, u_ptr, v_ptr)
730
731	24	if (showCallTree) call callTree_leave("btstep()")
732	24	call cpu_clock_end(id_clock_btstep)
733
734		! up = u + dt_pred*( u_bc_accel + u_accel_bt )
735	24	dt_pred = dt * CS%be
736	24	call cpu_clock_begin(id_clock_mom_update)
737
738	24	do concurrent (k=1:nz, J=Jsq:Jeq,i=is:ie)
739		vp(i,J,k) = G%mask2dCv(i,J) * (v_inst(i,J,k) + dt_pred * &
740	13354584	(v_bc_accel(i,J,k) + CS%v_accel_bt(i,J,k)))
741		enddo
742
743	24	do concurrent (k=1:nz, j=js:je, I=Isq:Ieq)
744		up(I,j,k) = G%mask2dCu(I,j) * (u_inst(I,j,k) + dt_pred * &
745	13245168	(u_bc_accel(I,j,k) + CS%u_accel_bt(I,j,k)))
746		enddo
747
748	24	call cpu_clock_end(id_clock_mom_update)
749
750	24	if (CS%debug) then
751		!$omp target update from(CS%CAu_pred, CS%CAv_pred, CS%PFu, CS%PFv)
752		!$omp target update from(CS%diffu, CS%diffv, CS%pbce)
753		!$omp target update from(CS%u_accel_bt, CS%v_accel_bt)
754		call MOM_accel_chksum("Predictor accel", CS%CAu_pred, CS%CAv_pred, CS%PFu, CS%PFv, &
755	0	CS%diffu, CS%diffv, G, GV, US, CS%pbce, CS%u_accel_bt, CS%v_accel_bt, symmetric=sym)
756
757		!$omp target update from(up, vp, h, uh, vh)
758	0	call uvchksum("Predictor 1 [uv]", up, vp, G%HI, haloshift=0, symmetric=sym, unscale=US%L_T_to_m_s)
759	0	call hchksum(h, "Predictor 1 h", G%HI, haloshift=1, unscale=GV%H_to_MKS)
760		call uvchksum("Predictor 1 [uv]h", uh, vh, G%HI,haloshift=2, &
761	0	symmetric=sym, unscale=GV%H_to_MKSUS%L_to_m2US%s_to_T)
762
763		! call MOM_state_chksum("Predictor 1", up, vp, h, uh, vh, G, GV, US, haloshift=1)
764
765		!$omp target update from(u_inst, v_inst)
766		call MOM_state_chksum("Predictor 1 init", u_inst, v_inst, h, uh, vh, G, GV, US, haloshift=1, &
767	0	symmetric=sym)
768	0	if (debug_redundant) then
769	0	call check_redundant("Predictor 1 up", up, vp, G, unscale=US%L_T_to_m_s)
770	0	call check_redundant("Predictor 1 uh", uh, vh, G, unscale=GV%H_to_MKSUS%L_to_m2US%s_to_T)
771		endif
772		endif
773
774		! up <- up + dt_pred d/dz visc d/dz up
775		! u_av <- u_av + dt_pred d/dz visc d/dz u_av
776	24	call cpu_clock_begin(id_clock_vertvisc)
777	24	if (CS%debug) then
778	0	call uvchksum("0 before vertvisc: [uv]p", up, vp, G%HI,haloshift=0, symmetric=sym, unscale=US%L_T_to_m_s)
779		endif
780
781	24	if (CS%fpmix) then
782		!$omp target update from(up, vp)
783	0	uold(:,:,:) = 0.0
784	0	vold(:,:,:) = 0.0
785	0	do k = 1, nz
786	0	do j = js , je
787	0	do I = Isq, Ieq
788	0	uold(I,j,k) = up(I,j,k)
789		enddo
790		enddo
791	0	do J = Jsq, Jeq
792	0	do i = is, ie
793	0	vold(i,J,k) = vp(i,J,k)
794		enddo
795		enddo
796		enddo
797		endif
798
799	24	call thickness_to_dz(h, tv, dz, G, GV, US, halo_size=1, do_offload=.true.)
800
801		call vertvisc_coef(up, vp, h, dz, forces, visc, tv, dt_pred, G, GV, US, CS%vertvisc_CSp, &
802	24	CS%OBC, VarMix)
803
804	24	if (CS%fpmix) then
805	0	hbl(:,:) = 0.0
806	0	if (ASSOCIATED(CS%KPP_CSp)) call KPP_get_BLD(CS%KPP_CSp, hbl, G, US, m_to_BLD_units=GV%m_to_H)
807	0	if (ASSOCIATED(CS%energetic_PBL_CSp)) &
808	0	call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, hbl, G, US, m_to_MLD_units=GV%m_to_H)
809
810		! lFPpost must be false in the predictor step to avoid averaging into the diagnostics
811	0	lFPpost = .false.
812		!$omp target update from(up, vp, h)
813		call vertFPmix(up, vp, uold, vold, hbl, h, forces, dt_pred, lFPpost, CS%Cemp_NL, &
814	0	G, GV, US, CS%vertvisc_CSp, CS%OBC, waves=waves)
815		!$omp target update to (up, vp)
816
817		call vertvisc(up, vp, h, forces, visc, dt_pred, CS%OBC, CS%AD_pred, CS%CDp, G, &
818	0	GV, US, CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot, fpmix=CS%fpmix, waves=waves)
819		else
820		call vertvisc(up, vp, h, forces, visc, dt_pred, CS%OBC, CS%AD_pred, CS%CDp, G, &
821	24	GV, US, CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot, waves=waves)
822
823		endif
824
825	24	if (showCallTree) call callTree_wayPoint("done with vertvisc (step_MOM_dyn_split_RK2)")
826	24	if (G%nonblocking_updates) then
827	0	call cpu_clock_end(id_clock_vertvisc)
828		!$omp target update from(up, vp)
829	0	call start_group_pass(CS%pass_uvp, G%Domain, clock=id_clock_pass)
830	0	call cpu_clock_begin(id_clock_vertvisc)
831		endif
832
833	24	if (CS%visc_rem_dt_bug) then
834	0	call vertvisc_remnant(visc, CS%visc_rem_u, CS%visc_rem_v, dt_pred, G, GV, US, CS%vertvisc_CSp)
835		else
836	24	call vertvisc_remnant(visc, CS%visc_rem_u, CS%visc_rem_v, dt, G, GV, US, CS%vertvisc_CSp)
837		endif
838
839	24	call cpu_clock_end(id_clock_vertvisc)
840
841	24	call do_group_pass(CS%pass_visc_rem, G%Domain, clock=id_clock_pass, omp_offload=.true.)
842
843	24	if (G%nonblocking_updates) then
844	0	call complete_group_pass(CS%pass_uvp, G%Domain, clock=id_clock_pass)
845		!$omp target update to(up, vp)
846		else
847	24	call do_group_pass(CS%pass_uvp, G%Domain, clock=id_clock_pass, omp_offload=.true.)
848		endif
849
850		! uh = u_av * h
851		! hp = h + dt * div . uh
852	24	call cpu_clock_begin(id_clock_continuity)
853		call continuity(up, vp, h, hp, uh, vh, dt, G, GV, US, CS%continuity_CSp, CS%OBC, pbv, &
854		uhbt=CS%uhbt, vhbt=CS%vhbt, visc_rem_u=CS%visc_rem_u, visc_rem_v=CS%visc_rem_v, &
855	24	u_cor=u_av, v_cor=v_av, BT_cont=CS%BT_cont)
856	24	call cpu_clock_end(id_clock_continuity)
857
858	24	if (showCallTree) call callTree_wayPoint("done with continuity (step_MOM_dyn_split_RK2)")
859
860	24	call do_group_pass(CS%pass_hp_uv, G%Domain, clock=id_clock_pass, omp_offload=.true.)
861
862	24	if (associated(CS%OBC)) then
863		!$omp target update from(u_av, v_av)
864
865	0	if (CS%debug) &
866	0	call uvchksum("Pre OBC avg [uv]", u_av, v_av, G%HI, haloshift=1, symmetric=sym, unscale=US%L_T_to_m_s)
867
868	0	call radiation_open_bdry_conds(CS%OBC, u_av, u_old_rad_OBC, v_av, v_old_rad_OBC, G, GV, US, dt_pred)
869
870	0	if (CS%debug) &
871	0	call uvchksum("Post OBC avg [uv]", u_av, v_av, G%HI, haloshift=1, symmetric=sym, unscale=US%L_T_to_m_s)
872
873		!$omp target update to(u_av, v_av)
874
875		! These should be done with a pass that excludes uh & vh.
876		! call do_group_pass(CS%pass_hp_uv, G%Domain, clock=id_clock_pass)
877	0	call pass_vector(u_av, v_av, G%Domain, halo=max(cor_stencil,vel_stencil), clock=id_clock_pass)
878		endif
879
880		! h_av = (h + hp)/2
881	3000	do concurrent (k=1:nz, j=js-cor_stencil:je+cor_stencil, i=is-cor_stencil:ie+cor_stencil)
882	14478264	h_av(i,j,k) = 0.5*(h(i,j,k) + hp(i,j,k))
883		enddo
884
885		! The correction phase of the time step starts here.
886	24	call enable_averages(dt, Time_local, CS%diag)
887
888		! Calculate a revised estimate of the free-surface height correction to be
889		! used in the next call to btstep. This call is at this point so that
890		! hp can be changed if CS%begw /= 0.
891		! eta_cor = ... (hidden inside CS%barotropic_CSp)
892	24	if (CS%BT_adj_corr_mass_src) then
893	24	call cpu_clock_begin(id_clock_btcalc)
894	24	call bt_mass_source(hp, eta_pred, .false., G, GV, CS%barotropic_CSp)
895	24	call cpu_clock_end(id_clock_btcalc)
896		endif
897
898	24	if (CS%begw /= 0.0) then
899		! hp <- (1-begw)h_in + begwhp
900		! Back up hp to the value it would have had after a time-step of
901		! begw*dt. hp is not used again until recalculated by continuity.
902	0	do concurrent (k=1:nz, j=js-2:je+2, i=is-2:ie+2)
903	0	hp(i,j,k) = (1.0-CS%begw)h(i,j,k) + CS%begwhp(i,j,k)
904		enddo
905
906		! PFu = d/dx M(hp,T,S)
907		! pbce = dM/deta
908	0	call cpu_clock_begin(id_clock_pres)
909		call PressureForce(hp, tv, CS%PFu, CS%PFv, G, GV, US, CS%PressureForce_CSp, &
910	0	CS%ALE_CSp, CS%ADp, p_surf, CS%pbce, CS%eta_PF)
911		! Stokes shear force contribution to pressure gradient
912	0	Use_Stokes_PGF = present(Waves)
913	0	if (Use_Stokes_PGF) then
914	0	Use_Stokes_PGF = associated(Waves)
915	0	if (Use_Stokes_PGF) Use_Stokes_PGF = Waves%Stokes_PGF
916	0	if (Use_Stokes_PGF) then
917		!$omp target update from(CS%PFu, CS%PFv, h)
918	0	call thickness_to_dz(h, tv, dz, G, GV, US, halo_size=1)
919	0	call Stokes_PGF(G, GV, US, dz, u_inst, v_inst, CS%PFu_Stokes, CS%PFv_Stokes, Waves)
920	0	if (.not.Waves%Passive_Stokes_PGF) then
921	0	do k=1,nz
922	0	do j=js,je ; do I=Isq,Ieq
923	0	CS%PFu(I,j,k) = CS%PFu(I,j,k) + CS%PFu_Stokes(I,j,k)
924		enddo ; enddo
925		enddo
926	0	do k=1,nz
927	0	do J=Jsq,Jeq ; do i=is,ie
928	0	CS%PFv(i,J,k) = CS%PFv(i,J,k) + CS%PFv_Stokes(i,J,k)
929		enddo ; enddo
930		enddo
931		endif
932		endif
933		endif
934	0	call cpu_clock_end(id_clock_pres)
935	0	if (showCallTree) call callTree_wayPoint("done with PressureForce[hp=(1-b).h+b.h] (step_MOM_dyn_split_RK2)")
936		endif
937
938		! TODO: Move p_surf handling outside of this function
939		!$omp target exit data map(delete: p_surf) if (associated(p_surf))
940
941	24	if (BT_cont_BT_thick) then
942		call btcalc(h, G, GV, CS%barotropic_CSp, CS%BT_cont%h_u, CS%BT_cont%h_v, &
943	24	OBC=CS%OBC)
944	24	if (showCallTree) call callTree_wayPoint("done with btcalc[BT_cont_BT_thick] (step_MOM_dyn_split_RK2)")
945		endif
946
947	24	if (CS%debug) then
948		!$omp target update from(up, vp, hp, uh, vh)
949	0	call MOM_state_chksum("Predictor ", up, vp, hp, uh, vh, G, GV, US, symmetric=sym)
950		!$omp target update from(u_av, v_av, h_av)
951	0	call uvchksum("Predictor avg [uv]", u_av, v_av, G%HI, haloshift=1, symmetric=sym, unscale=US%L_T_to_m_s)
952	0	call hchksum(h_av, "Predictor avg h", G%HI, haloshift=2, unscale=GV%H_to_MKS)
953		! call MOM_state_chksum("Predictor avg ", u_av, v_av, h_av, uh, vh, G, GV, US)
954	0	if (debug_redundant) then
955	0	call check_redundant("Predictor up ", up, vp, G, unscale=US%L_T_to_m_s)
956	0	call check_redundant("Predictor uh ", uh, vh, G, unscale=GV%H_to_MKSUS%L_to_m2US%s_to_T)
957		endif
958		endif
959
960		! diffu = horizontal viscosity terms (u_av)
961	24	call cpu_clock_begin(id_clock_horvisc)
962		call horizontal_viscosity(u_av, v_av, h_av, uh, vh, CS%diffu, CS%diffv, &
963		MEKE, Varmix, G, GV, US, CS%hor_visc, hor_visc_nkblock(CS%hor_visc, GV), tv, dt, &
964		OBC=CS%OBC, BT=CS%barotropic_CSp, TD=thickness_diffuse_CSp, &
965	24	ADp=CS%ADp, hu_cont=CS%BT_cont%h_u, hv_cont=CS%BT_cont%h_v, STOCH=STOCH)
966
967	24	call cpu_clock_end(id_clock_horvisc)
968	24	if (showCallTree) call callTree_wayPoint("done with horizontal_viscosity (step_MOM_dyn_split_RK2)")
969
970		! CAu = -(f+zeta_av)/h_av vh + d/dx KE_av
971	24	call cpu_clock_begin(id_clock_Cor)
972		call CorAdCalc(u_av, v_av, h_av, uh, vh, CS%CAu, CS%CAv, CS%OBC, CS%ADp, &
973	24	G, GV, US, CS%CoriolisAdv, pbv, Waves=Waves)
974	24	call cpu_clock_end(id_clock_Cor)
975
976	24	if (showCallTree) call callTree_wayPoint("done with CorAdCalc (step_MOM_dyn_split_RK2)")
977
978		! Calculate the momentum forcing terms for the barotropic equations.
979
980		! u_bc_accel = CAu + PFu + diffu(u[n-1])
981	24	call cpu_clock_begin(id_clock_btforce)
982
983	24	do concurrent (k=1:nz, j=js:je, I=Isq:Ieq)
984	13245168	u_bc_accel(I,j,k) = (CS%Cau(I,j,k) + CS%PFu(I,j,k)) + CS%diffu(I,j,k)
985		enddo
986
987	24	do concurrent (k=1:nz, J=Jsq:Jeq,i=is:ie)
988	13354584	v_bc_accel(i,J,k) = (CS%Cav(i,J,k) + CS%PFv(i,J,k)) + CS%diffv(i,J,k)
989		enddo
990
991	24	if (associated(CS%OBC)) then
992		!$omp target update from(u_bc_accel, v_bc_accel)
993	0	call open_boundary_zero_normal_flow(CS%OBC, G, GV, u_bc_accel, v_bc_accel)
994		!$omp target update to(u_bc_accel, v_bc_accel)
995		endif
996	24	call cpu_clock_end(id_clock_btforce)
997
998	24	if (CS%debug) then
999		!$omp target update from(CS%CAu, CS%CAv)
1000		!$omp target update from(CS%PFu, CS%PFv, CS%pbce)
1001		!$omp target update from(CS%diffu, CS%diffv)
1002		!$omp target update from(u_bc_accel, v_bc_accel)
1003		call MOM_accel_chksum("corr pre-btstep accel", CS%CAu, CS%CAv, CS%PFu, CS%PFv, &
1004		CS%diffu, CS%diffv, G, GV, US, CS%pbce, u_bc_accel, v_bc_accel, &
1005	0	symmetric=sym)
1006	0	if (debug_redundant) then
1007	0	call check_redundant("corr pre-btstep CS%CA ", CS%CAu, CS%CAv, G, unscale=US%L_T2_to_m_s2)
1008	0	call check_redundant("corr pre-btstep CS%PF ", CS%PFu, CS%PFv, G, unscale=US%L_T2_to_m_s2)
1009	0	call check_redundant("corr pre-btstep CS%diff ", CS%diffu, CS%diffv, G, unscale=US%L_T2_to_m_s2)
1010	0	call check_redundant("corr pre-btstep u_bc_accel ", u_bc_accel, v_bc_accel, G, unscale=US%L_T2_to_m_s2)
1011		endif
1012		endif
1013
1014		! u_accel_bt = layer accelerations due to barotropic solver
1015		! pbce = dM/deta
1016	24	call cpu_clock_begin(id_clock_btstep)
1017	24	if (CS%BT_use_layer_fluxes) then
1018	24	uh_ptr => uh ; vh_ptr => vh ; u_ptr => u_av ; v_ptr => v_av
1019		endif
1020
1021	24	if (showCallTree) call callTree_enter("btstep(), MOM_barotropic.F90")
1022		! This is the corrector step call to btstep.
1023		call btstep(u_inst, v_inst, eta, dt, u_bc_accel, v_bc_accel, forces, CS%pbce, CS%eta_PF, u_av, v_av, &
1024		CS%u_accel_bt, CS%v_accel_bt, eta_pred, CS%uhbt, CS%vhbt, G, GV, US, &
1025		CS%barotropic_CSp, CS%visc_rem_u, CS%visc_rem_v, SpV_avg, CS%ADp, CS%OBC, CS%BT_cont, &
1026	24	eta_PF_start, taux_bot, tauy_bot, uh_ptr, vh_ptr, u_ptr, v_ptr, etaav=eta_av)
1027
1028	24	if (CS%id_deta_dt>0) then
1029		!$omp target update from(eta, eta_pred)
1030	0	do j=js,je ; do i=is,ie ; deta_dt(i,j) = (eta_pred(i,j) - eta(i,j))*Idt_bc ; enddo ; enddo
1031		endif
1032
1033	24	do concurrent (j=js:je, i=is:ie)
1034	175704	eta(i,j) = eta_pred(i,j)
1035		enddo
1036
1037	24	call cpu_clock_end(id_clock_btstep)
1038	24	if (showCallTree) call callTree_leave("btstep()")
1039
1040	24	if (CS%debug .and. debug_redundant) then
1041		!$omp target update from(CS%u_accel_bt, CS%v_accel_bt)
1042	0	call check_redundant("u_accel_bt ", CS%u_accel_bt, CS%v_accel_bt, G, unscale=US%L_T2_to_m_s2)
1043		endif
1044
1045		! u = u + dt*( u_bc_accel + u_accel_bt )
1046	24	call cpu_clock_begin(id_clock_mom_update)
1047	24	do concurrent (k=1:nz, j=js:je, I=Isq:Ieq)
1048		u_inst(I,j,k) = G%mask2dCu(I,j) * (u_inst(I,j,k) + dt * &
1049	13245168	(u_bc_accel(I,j,k) + CS%u_accel_bt(I,j,k)))
1050		enddo
1051	24	do concurrent (k=1:nz, J=Jsq:Jeq, i=is:ie)
1052		v_inst(i,J,k) = G%mask2dCv(i,J) * (v_inst(i,J,k) + dt * &
1053	13354584	(v_bc_accel(i,J,k) + CS%v_accel_bt(i,J,k)))
1054		enddo
1055	24	call cpu_clock_end(id_clock_mom_update)
1056
1057	24	if (CS%debug) then
1058		!$omp target update from(CS%CAu, CS%CAv)
1059		!$omp target update from(CS%PFu, CS%PFv, CS%pbce)
1060		!$omp target update from(CS%diffu, CS%diffv)
1061		!$omp target update from(CS%u_accel_bt, CS%v_accel_bt)
1062		call MOM_accel_chksum("Corrector accel", CS%CAu, CS%CAv, CS%PFu, CS%PFv, &
1063		CS%diffu, CS%diffv, G, GV, US, CS%pbce, CS%u_accel_bt, CS%v_accel_bt, &
1064	0	symmetric=sym)
1065
1066		!$omp target update from(u_inst, v_inst, h, uh, vh)
1067	0	call uvchksum("Corrector 1 [uv]", u_inst, v_inst, G%HI, haloshift=0, symmetric=sym, unscale=US%L_T_to_m_s)
1068	0	call hchksum(h, "Corrector 1 h", G%HI, haloshift=1, unscale=GV%H_to_MKS)
1069		call uvchksum("Corrector 1 [uv]h", uh, vh, G%HI, haloshift=2, &
1070	0	symmetric=sym, unscale=GV%H_to_MKSUS%L_to_m2US%s_to_T)
1071
1072		! call MOM_state_chksum("Corrector 1", u_inst, v_inst, h, uh, vh, G, GV, US, haloshift=1)
1073		endif
1074
1075		! u <- u + dt d/dz visc d/dz u
1076		! u_av <- u_av + dt d/dz visc d/dz u_av
1077	24	call cpu_clock_begin(id_clock_vertvisc)
1078
1079	24	if (CS%fpmix) then
1080		!$omp target update from(u_inst, v_inst)
1081	0	uold(:,:,:) = 0.0
1082	0	vold(:,:,:) = 0.0
1083	0	do k = 1, nz
1084	0	do j = js , je
1085	0	do I = Isq, Ieq
1086	0	uold(I,j,k) = u_inst(I,j,k)
1087		enddo
1088		enddo
1089	0	do J = Jsq, Jeq
1090	0	do i = is, ie
1091	0	vold(i,J,k) = v_inst(i,J,k)
1092		enddo
1093		enddo
1094		enddo
1095		endif
1096
1097	24	call thickness_to_dz(h, tv, dz, G, GV, US, halo_size=1, do_offload=.true.)
1098
1099	24	call vertvisc_coef(u_inst, v_inst, h, dz, forces, visc, tv, dt, G, GV, US, CS%vertvisc_CSp, CS%OBC, VarMix)
1100
1101	24	if (CS%fpmix) then
1102	0	lFPpost = .true.
1103		!$omp target update from(u_inst, v_inst, h)
1104		call vertFPmix(u_inst, v_inst, uold, vold, hbl, h, forces, dt, lFPpost, CS%Cemp_NL, &
1105	0	G, GV, US, CS%vertvisc_CSp, CS%OBC, Waves=Waves)
1106		!$omp target update to(u_inst, v_inst)
1107
1108		call vertvisc(u_inst, v_inst, h, forces, visc, dt, CS%OBC, CS%ADp, CS%CDp, G, GV, US, &
1109	0	CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot, fpmix=CS%fpmix, waves=waves)
1110		else
1111		call vertvisc(u_inst, v_inst, h, forces, visc, dt, CS%OBC, CS%ADp, CS%CDp, G, GV, US, &
1112	24	CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot, waves=waves)
1113		endif
1114
1115	24	if (G%nonblocking_updates) then
1116		!$omp target update from(u_inst, v_inst)
1117	0	call cpu_clock_end(id_clock_vertvisc)
1118	0	call start_group_pass(CS%pass_uv, G%Domain, clock=id_clock_pass)
1119	0	call cpu_clock_begin(id_clock_vertvisc)
1120		endif
1121	24	call vertvisc_remnant(visc, CS%visc_rem_u, CS%visc_rem_v, dt, G, GV, US, CS%vertvisc_CSp)
1122
1123	24	call cpu_clock_end(id_clock_vertvisc)
1124	24	if (showCallTree) call callTree_wayPoint("done with vertvisc (step_MOM_dyn_split_RK2)")
1125
1126		! Later, h_av = (h_in + h_out)/2, but for now use h_av to store h_in.
1127	3000	do concurrent (k=1:nz, j=js-cor_stencil:je+cor_stencil, i=is-cor_stencil:ie+cor_stencil)
1128	14478264	h_av(i,j,k) = h(i,j,k)
1129		enddo
1130
1131	24	call do_group_pass(CS%pass_visc_rem, G%Domain, clock=id_clock_pass, omp_offload=.true.)
1132
1133	24	if (G%nonblocking_updates) then
1134	0	call complete_group_pass(CS%pass_uv, G%Domain, clock=id_clock_pass)
1135		!$omp target update to(u_inst, v_inst)
1136		else
1137	24	call do_group_pass(CS%pass_uv, G%Domain, clock=id_clock_pass, omp_offload=.true.)
1138		endif
1139
1140		! uh = u_av * h
1141		! h = h + dt * div . uh
1142		! u_av and v_av adjusted so their mass transports match uhbt and vhbt.
1143	24	call cpu_clock_begin(id_clock_continuity)
1144	3096	do concurrent (k=1:nz, j=G%jsd:G%jed, i=G%isd:G%ied)
1145	15879192	h_tmp(i,j,k) = h(i,j,k)
1146		enddo
1147
1148		call continuity(u_inst, v_inst, h_tmp, h, uh, vh, dt, G, GV, US, CS%continuity_CSp, CS%OBC, pbv, &
1149		uhbt=CS%uhbt, vhbt=CS%vhbt, visc_rem_u=CS%visc_rem_u, visc_rem_v=CS%visc_rem_v, &
1150	24	u_cor=u_av, v_cor=v_av)
1151	24	call cpu_clock_end(id_clock_continuity)
1152
1153	24	call do_group_pass(CS%pass_h, G%Domain, clock=id_clock_pass, omp_offload=.true.)
1154
1155
1156		! Whenever thickness changes let the diag manager know, target grids
1157		! for vertical remapping may need to be regenerated.
1158	24	call diag_update_remap_grids(CS%diag)
1159	24	if (showCallTree) call callTree_wayPoint("done with continuity (step_MOM_dyn_split_RK2)")
1160
1161	24	if (G%nonblocking_updates) then
1162		!$omp target update from(u_av, v_av, uh, vh)
1163	0	call start_group_pass(CS%pass_av_uvh, G%Domain, clock=id_clock_pass)
1164		else
1165	24	call do_group_pass(CS%pass_av_uvh, G%domain, clock=id_clock_pass, omp_offload=.true.)
1166		endif
1167
1168	24	if (associated(CS%OBC)) then
1169		!### I suspect that there is a bug here when u_inst is compared with a previous value of u_av
1170		! to estimate the dominant outward group velocity, but a fix is not available yet.
1171		!$omp target update from(u_inst, v_inst)
1172	0	call radiation_open_bdry_conds(CS%OBC, u_inst, u_old_rad_OBC, v_inst, v_old_rad_OBC, G, GV, US, dt)
1173		endif
1174
1175		! h_av = (h_in + h_out)/2 . Going in to this line, h_av = h_in.
1176	3000	do concurrent (k=1:nz, j=js-cor_stencil:je+cor_stencil, i=is-cor_stencil:ie+cor_stencil)
1177	14478264	h_av(i,j,k) = 0.5*(h_av(i,j,k) + h(i,j,k))
1178		enddo
1179
1180	24	if (G%nonblocking_updates) then
1181	0	call complete_group_pass(CS%pass_av_uvh, G%Domain, clock=id_clock_pass)
1182		!$omp target update to(u_av, v_av, uh, vh)
1183		endif
1184
1185	3024	do concurrent (k=1:nz, j=js-cor_stencil:je+cor_stencil, I=Isq-cor_stencil:Ieq+cor_stencil)
1186	14595024	uhtr(I,j,k) = uhtr(I,j,k) + uh(I,j,k) * dt
1187		enddo
1188
1189	3000	do concurrent (k=1:nz, J=Jsq-cor_stencil:Jeq+cor_stencil, i=is-cor_stencil:ie+cor_stencil)
1190	14704440	vhtr(i,J,k) = vhtr(i,J,k) + vh(i,J,k) * dt
1191		enddo
1192
1193		!$omp target exit data map(release: u_bc_accel, v_bc_accel, eta_pred, uh_in, vh_in)
1194		!$omp target exit data map(delete: hp, up, vp, dz, h_tmp)
1195
1196	24	if (associated(CS%OBC)) then
1197	0	call update_segment_thickness_reservoirs(G, GV, uhtr, vhtr, h, CS%OBC)
1198		endif
1199
1200	24	if (CS%store_CAu) then
1201		! Calculate a predictor-step estimate of the Coriolis and momentum advection terms
1202		! for use in the next time step, possibly after it has been vertically remapped.
1203	24	call cpu_clock_begin(id_clock_Cor)
1204	24	call disable_averaging(CS%diag) ! These calculations should not be used for diagnostics.
1205		! CAu = -(f+zeta_av)/h_av vh + d/dx KE_av
1206		call CorAdCalc(u_av, v_av, h_av, uh, vh, CS%CAu_pred, CS%CAv_pred, CS%OBC, CS%AD_pred, &
1207	24	G, GV, US, CS%CoriolisAdv, pbv, Waves=Waves)
1208	24	CS%CAu_pred_stored = .true.
1209	24	call enable_averages(dt, Time_local, CS%diag) ! Reenable the averaging
1210	24	call cpu_clock_end(id_clock_Cor)
1211	24	if (showCallTree) call callTree_wayPoint("done with CorAdCalc (step_MOM_dyn_split_RK2)")
1212		else
1213	0	CS%CAu_pred_stored = .false.
1214		endif
1215
1216		! The time-averaged free surface height has already been set by the last call to btstep.
1217
1218		! Deallocate this memory to avoid a memory leak. ### We should revisit how this array is declared. -RWH
1219	24	if (dyn_p_surf .and. associated(eta_PF_start)) deallocate(eta_PF_start)
1220
1221		! Here various terms used in to update the momentum equations are
1222		! offered for time averaging.
1223	24	if (CS%id_PFu > 0) call post_data(CS%id_PFu, CS%PFu, CS%diag)
1224	24	if (CS%id_PFv > 0) call post_data(CS%id_PFv, CS%PFv, CS%diag)
1225	24	if (CS%id_CAu > 0) call post_data(CS%id_CAu, CS%CAu, CS%diag)
1226	24	if (CS%id_CAv > 0) call post_data(CS%id_CAv, CS%CAv, CS%diag)
1227
1228		! Here the thickness fluxes are offered for time averaging.
1229	24	if (CS%id_uh > 0) call post_data(CS%id_uh, uh, CS%diag)
1230	24	if (CS%id_vh > 0) call post_data(CS%id_vh, vh, CS%diag)
1231	24	if (CS%id_uav > 0) call post_data(CS%id_uav, u_av, CS%diag)
1232	24	if (CS%id_vav > 0) call post_data(CS%id_vav, v_av, CS%diag)
1233	24	if (CS%id_u_BT_accel > 0) call post_data(CS%id_u_BT_accel, CS%u_accel_bt, CS%diag)
1234	24	if (CS%id_v_BT_accel > 0) call post_data(CS%id_v_BT_accel, CS%v_accel_bt, CS%diag)
1235
1236		! Calculate effective areas and post data
1237	24	if (CS%id_ueffA > 0) then
1238	0	ueffA(:,:,:) = 0
1239	0	do k=1,nz ; do j=js,je ; do I=Isq,Ieq
1240	0	if (abs(up(I,j,k)) > 0.) ueffA(I,j,k) = uh(I,j,k) / up(I,j,k)
1241		enddo ; enddo ; enddo
1242	0	call post_data(CS%id_ueffA, ueffA, CS%diag)
1243		endif
1244
1245	24	if (CS%id_veffA > 0) then
1246	0	veffA(:,:,:) = 0
1247	0	do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
1248	0	if (abs(vp(i,J,k)) > 0.) veffA(i,J,k) = vh(i,J,k) / vp(i,J,k)
1249		enddo ; enddo ; enddo
1250	0	call post_data(CS%id_veffA, veffA, CS%diag)
1251		endif
1252
1253		!$omp target exit data map(delete: up, vp)
1254
1255		! Diagnostics of the fractional thicknesses times momentum budget terms
1256		! 3D diagnostics hf_PFu etc. are commented because there is no clarity on proper remapping grid option.
1257		! The code is retained for debugging purposes in the future.
1258		!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)
1259		!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)
1260		!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)
1261		!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)
1262		!if (CS%id_hf_u_BT_accel > 0) &
1263		! call post_product_u(CS%id_hf_u_BT_accel, CS%u_accel_bt, CS%ADp%diag_hfrac_u, G, nz, CS%diag)
1264		!if (CS%id_hf_v_BT_accel > 0) &
1265		! call post_product_v(CS%id_hf_v_BT_accel, CS%v_accel_bt, CS%ADp%diag_hfrac_v, G, nz, CS%diag)
1266
1267		! Diagnostics for the vertical sum of layer thickness x prssure force accelerations
1268	24	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)
1269	24	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)
1270
1271		! Diagnostics for thickness-weighted vertically averaged prssure force accelerations
1272	24	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)
1273	24	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)
1274
1275		! Diagnostics for thickness x prssure force accelerations
1276	24	if (CS%id_h_PFu > 0) call post_product_u(CS%id_h_PFu, CS%PFu, CS%ADp%diag_hu, G, nz, CS%diag)
1277	24	if (CS%id_h_PFv > 0) call post_product_v(CS%id_h_PFv, CS%PFv, CS%ADp%diag_hv, G, nz, CS%diag)
1278
1279		! Diagnostics of Coriolis acceleratations
1280	24	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)
1281	24	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)
1282	24	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)
1283	24	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)
1284	24	if (CS%id_h_CAu > 0) call post_product_u(CS%id_h_CAu, CS%CAu, CS%ADp%diag_hu, G, nz, CS%diag)
1285	24	if (CS%id_h_CAv > 0) call post_product_v(CS%id_h_CAv, CS%CAv, CS%ADp%diag_hv, G, nz, CS%diag)
1286
1287		! Diagnostics of barotropic solver acceleratations
1288	24	if (CS%id_intz_u_BT_accel_2d > 0) &
1289	0	call post_product_sum_u(CS%id_intz_u_BT_accel_2d, CS%u_accel_bt, CS%ADp%diag_hu, G, nz, CS%diag)
1290	24	if (CS%id_intz_v_BT_accel_2d > 0) &
1291	0	call post_product_sum_v(CS%id_intz_v_BT_accel_2d, CS%v_accel_bt, CS%ADp%diag_hv, G, nz, CS%diag)
1292	24	if (CS%id_hf_u_BT_accel_2d > 0) &
1293	0	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)
1294	24	if (CS%id_hf_v_BT_accel_2d > 0) &
1295	0	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)
1296	24	if (CS%id_h_u_BT_accel > 0) &
1297	0	call post_product_u(CS%id_h_u_BT_accel, CS%u_accel_bt, CS%ADp%diag_hu, G, nz, CS%diag)
1298	24	if (CS%id_h_v_BT_accel > 0) &
1299	0	call post_product_v(CS%id_h_v_BT_accel, CS%v_accel_bt, CS%ADp%diag_hv, G, nz, CS%diag)
1300
1301		! Diagnostics for momentum budget terms multiplied by visc_rem_[uv],
1302	24	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)
1303	24	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)
1304	24	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)
1305	24	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)
1306	24	if (CS%id_u_BT_accel_visc_rem > 0) &
1307	0	call post_product_u(CS%id_u_BT_accel_visc_rem, CS%u_accel_bt, CS%ADp%visc_rem_u, G, nz, CS%diag)
1308	24	if (CS%id_v_BT_accel_visc_rem > 0) &
1309	0	call post_product_v(CS%id_v_BT_accel_visc_rem, CS%v_accel_bt, CS%ADp%visc_rem_v, G, nz, CS%diag)
1310
1311		! Diagnostics related to changes in eta
1312	24	if (CS%id_deta_dt > 0) call post_data(CS%id_deta_dt, deta_dt, CS%diag)
1313
1314	24	if (CS%debug) then
1315		!$omp target update from(u_inst, v_inst, h, uh, vh, u_av, v_av, h_av)
1316	0	call MOM_state_chksum("Corrector ", u_inst, v_inst, h, uh, vh, G, GV, US, symmetric=sym)
1317	0	call uvchksum("Corrector avg [uv]", u_av, v_av, G%HI, haloshift=1, symmetric=sym, unscale=US%L_T_to_m_s)
1318	0	call hchksum(h_av, "Corrector avg h", G%HI, haloshift=1, unscale=GV%H_to_MKS)
1319		! call MOM_state_chksum("Corrector avg ", u_av, v_av, h_av, uh, vh, G, GV, US)
1320		endif
1321
1322	24	if (showCallTree) call callTree_leave("step_MOM_dyn_split_RK2()")
1323
1324	24	end subroutine step_MOM_dyn_split_RK2
1325
1326		!> This subroutine sets up any auxiliary restart variables that are specific
1327		!! to the split-explicit time stepping scheme. All variables registered here should
1328		!! have the ability to be recreated if they are not present in a restart file.
1329	1	subroutine register_restarts_dyn_split_RK2(HI, GV, US, param_file, CS, restart_CS, uh, vh)
1330		type(hor_index_type), intent(in) :: HI !< Horizontal index structure
1331		type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
1332		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1333		type(param_file_type), intent(in) :: param_file !< parameter file
1334		type(MOM_dyn_split_RK2_CS), intent(inout) :: CS !< module control structure
1335		type(MOM_restart_CS), intent(inout) :: restart_CS !< MOM restart control structure
1336		real, dimension(SZIB_(HI),SZJ_(HI),SZK_(GV)), &
1337		target, intent(inout) :: uh !< zonal volume or mass transport [H L2 T-1 ~> m3 s-1 or kg s-1]
1338		real, dimension(SZI_(HI),SZJB_(HI),SZK_(GV)), &
1339		target, intent(inout) :: vh !< merid volume or mass transport [H L2 T-1 ~> m3 s-1 or kg s-1]
1340
1341		character(len=40) :: mdl = "MOM_dynamics_split_RK2" ! This module's name.
1342	41	type(vardesc) :: vd(2)
1343		character(len=48) :: thickness_units, flux_units
1344
1345		integer :: isd, ied, jsd, jed, nz, IsdB, IedB, JsdB, JedB
1346
1347	1	isd = HI%isd ; ied = HI%ied ; jsd = HI%jsd ; jed = HI%jed ; nz = GV%ke
1348	1	IsdB = HI%IsdB ; IedB = HI%IedB ; JsdB = HI%JsdB ; JedB = HI%JedB
1349
1350		! TODO: Are these initializations necessary? If not, then we can do
1351		! map(alloc:) rather than map(to:)
1352	663076	ALLOC_(CS%diffu(IsdB:IedB,jsd:jed,nz)) ; CS%diffu(:,:,:) = 0.0
1353	667651	ALLOC_(CS%diffv(isd:ied,JsdB:JedB,nz)) ; CS%diffv(:,:,:) = 0.0
1354		!$omp target enter data map(to: CS%diffu, CS%diffv)
1355	663076	ALLOC_(CS%CAu(IsdB:IedB,jsd:jed,nz)) ; CS%CAu(:,:,:) = 0.0
1356	667651	ALLOC_(CS%CAv(isd:ied,JsdB:JedB,nz)) ; CS%CAv(:,:,:) = 0.0
1357		!$omp target enter data map(to: CS%CAu, CS%CAv)
1358	663076	ALLOC_(CS%CAu_pred(IsdB:IedB,jsd:jed,nz)) ; CS%CAu_pred(:,:,:) = 0.0
1359	667651	ALLOC_(CS%CAv_pred(isd:ied,JsdB:JedB,nz)) ; CS%CAv_pred(:,:,:) = 0.0
1360		!$omp target enter data map(to: CS%CAu_pred, CS%CAv_pred)
1361	663076	ALLOC_(CS%PFu(IsdB:IedB,jsd:jed,nz)) ; CS%PFu(:,:,:) = 0.0
1362	667651	ALLOC_(CS%PFv(isd:ied,JsdB:JedB,nz)) ; CS%PFv(:,:,:) = 0.0
1363		!$omp target enter data map(to: CS%PFu, CS%PFv)
1364	8773	ALLOC_(CS%eta(isd:ied,jsd:jed)) ; CS%eta(:,:) = 0.0
1365	663076	ALLOC_(CS%u_av(IsdB:IedB,jsd:jed,nz)) ; CS%u_av(:,:,:) = 0.0
1366	667651	ALLOC_(CS%v_av(isd:ied,JsdB:JedB,nz)) ; CS%v_av(:,:,:) = 0.0
1367	657976	ALLOC_(CS%h_av(isd:ied,jsd:jed,nz)) ; CS%h_av(:,:,:) = GV%Angstrom_H
1368		!$omp target enter data map(to: CS%eta, CS%u_av, CS%v_av, CS%h_av)
1369
1370	1	thickness_units = get_thickness_units(GV)
1371	1	flux_units = get_flux_units(GV)
1372
1373		call get_param(param_file, mdl, "STORE_CORIOLIS_ACCEL", CS%store_CAu, &
1374		"If true, calculate the Coriolis accelerations at the end of each "//&
1375		"timestep for use in the predictor step of the next split RK2 timestep.", &
1376	1	default=.true., do_not_log=.true.)
1377
1378	1	if (GV%Boussinesq) then
1379		call register_restart_field(CS%eta, "sfc", .false., restart_CS, &
1380	1	longname="Free surface Height", units=thickness_units, conversion=GV%H_to_mks)
1381		else
1382		call register_restart_field(CS%eta, "p_bot", .false., restart_CS, &
1383	0	longname="Bottom Pressure", units=thickness_units, conversion=GV%H_to_mks)
1384		endif
1385
1386		! These are needed, either to calculate CAu and CAv or to calculate the velocity anomalies in
1387		! the barotropic solver's Coriolis terms.
1388	6	vd(1) = var_desc("u2", "m s-1", "Auxiliary Zonal velocity", 'u', 'L')
1389	6	vd(2) = var_desc("v2", "m s-1", "Auxiliary Meridional velocity", 'v', 'L')
1390		call register_restart_pair(CS%u_av, CS%v_av, vd(1), vd(2), .false., restart_CS, &
1391	1	conversion=US%L_T_to_m_s)
1392
1393	1	if (CS%store_CAu) then
1394	6	vd(1) = var_desc("CAu", "m s-2", "Zonal Coriolis and advactive acceleration", 'u', 'L')
1395	6	vd(2) = var_desc("CAv", "m s-2", "Meridional Coriolis and advactive acceleration", 'v', 'L')
1396		call register_restart_pair(CS%CAu_pred, CS%CAv_pred, vd(1), vd(2), .false., restart_CS, &
1397	1	conversion=US%L_T2_to_m_s2)
1398		else
1399		call register_restart_field(CS%h_av, "h2", .false., restart_CS, &
1400	0	longname="Auxiliary Layer Thickness", units=thickness_units, conversion=GV%H_to_mks)
1401
1402	0	vd(1) = var_desc("uh", flux_units, "Zonal thickness flux", 'u', 'L')
1403	0	vd(2) = var_desc("vh", flux_units, "Meridional thickness flux", 'v', 'L')
1404		call register_restart_pair(uh, vh, vd(1), vd(2), .false., restart_CS, &
1405	0	conversion=GV%H_to_MKSUS%L_to_m2US%s_to_T)
1406		endif
1407
1408	6	vd(1) = var_desc("diffu", "m s-2", "Zonal horizontal viscous acceleration", 'u', 'L')
1409	6	vd(2) = var_desc("diffv", "m s-2", "Meridional horizontal viscous acceleration", 'v', 'L')
1410		call register_restart_pair(CS%diffu, CS%diffv, vd(1), vd(2), .false., restart_CS, &
1411	1	conversion=US%L_T2_to_m_s2)
1412
1413	1	call register_barotropic_restarts(HI, GV, US, param_file, CS%barotropic_CSp, restart_CS)
1414
1415	22	end subroutine register_restarts_dyn_split_RK2
1416
1417		!> This subroutine does remapping for the auxiliary restart variables that are used
1418		!! with the split RK2 time stepping scheme.
1419	0	subroutine remap_dyn_split_RK2_aux_vars(G, GV, CS, h_old_u, h_old_v, h_new_u, h_new_v, ALE_CSp)
1420		type(ocean_grid_type), intent(inout) :: G !< ocean grid structure
1421		type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
1422		type(MOM_dyn_split_RK2_CS), pointer :: CS !< module control structure
1423		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
1424		intent(in) :: h_old_u !< Source grid thickness at zonal
1425		!! velocity points [H ~> m or kg m-2]
1426		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
1427		intent(in) :: h_old_v !< Source grid thickness at meridional
1428		!! velocity points [H ~> m or kg m-2]
1429		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
1430		intent(in) :: h_new_u !< Destination grid thickness at zonal
1431		!! velocity points [H ~> m or kg m-2]
1432		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
1433		intent(in) :: h_new_v !< Destination grid thickness at meridional
1434		!! velocity points [H ~> m or kg m-2]
1435		type(ALE_CS), pointer :: ALE_CSp !< ALE control structure to use when remapping
1436
1437	0	if (.not.CS%remap_aux) return
1438
1439	0	if (CS%store_CAu) then
1440	0	call ALE_remap_velocities(ALE_CSp, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, CS%u_av, CS%v_av)
1441	0	call pass_vector(CS%u_av, CS%v_av, G%Domain, complete=.false.)
1442	0	call ALE_remap_velocities(ALE_CSp, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, CS%CAu_pred, CS%CAv_pred)
1443	0	call pass_vector(CS%CAu_pred, CS%CAv_pred, G%Domain, complete=.true.)
1444		endif
1445
1446	0	call ALE_remap_velocities(ALE_CSp, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, CS%diffu, CS%diffv)
1447
1448		end subroutine remap_dyn_split_RK2_aux_vars
1449
1450		!> Initializes aspects of the dyn_split_RK2 that depend on diabatic processes.
1451		!! Needed when BLDs are used in the dynamics.
1452	0	subroutine init_dyn_split_RK2_diabatic(diabatic_CSp, CS)
1453		type(diabatic_CS), intent(in) :: diabatic_CSp !< diabatic structure
1454		type(MOM_dyn_split_RK2_CS), pointer :: CS !< module control structure
1455
1456	0	call extract_diabatic_member(diabatic_CSp, KPP_CSp=CS%KPP_CSp)
1457	0	call extract_diabatic_member(diabatic_CSp, energetic_PBL_CSp=CS%energetic_PBL_CSp)
1458
1459	0	end subroutine init_dyn_split_RK2_diabatic
1460
1461		!> This subroutine initializes all of the variables that are used by this
1462		!! dynamic core, including diagnostics and the cpu clocks.
1463	1	subroutine initialize_dyn_split_RK2(u, v, h, tv, uh, vh, eta, Time, G, GV, US, param_file, &
1464		diag, CS, HA_CSp, restart_CS, dt, Accel_diag, Cont_diag, MIS, &
1465		VarMix, MEKE, thickness_diffuse_CSp, &
1466		OBC, update_OBC_CSp, ALE_CSp, set_visc, &
1467		visc, dirs, ntrunc, pbv, calc_dtbt, cont_stencil, dyn_h_stencil)
1468		type(ocean_grid_type), intent(inout) :: G !< ocean grid structure
1469		type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
1470		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1471		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
1472		intent(inout) :: u !< zonal velocity [L T-1 ~> m s-1]
1473		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
1474		intent(inout) :: v !< merid velocity [L T-1 ~> m s-1]
1475		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1476		intent(inout) :: h !< layer thickness [H ~> m or kg m-2]
1477		type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic type
1478		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
1479		target, intent(inout) :: uh !< zonal volume/mass transport [H L2 T-1 ~> m3 s-1 or kg s-1]
1480		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
1481		target, intent(inout) :: vh !< merid volume/mass transport [H L2 T-1 ~> m3 s-1 or kg s-1]
1482		real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: eta !< free surface height or column mass [H ~> m or kg m-2]
1483		type(time_type), target, intent(in) :: Time !< current model time
1484		type(param_file_type), intent(in) :: param_file !< parameter file for parsing
1485		type(diag_ctrl), target, intent(inout) :: diag !< to control diagnostics
1486		type(MOM_dyn_split_RK2_CS), pointer :: CS !< module control structure
1487		type(harmonic_analysis_CS), pointer :: HA_CSp !< A pointer to the control structure of the
1488		!! harmonic analysis module
1489		type(MOM_restart_CS), intent(inout) :: restart_CS !< MOM restart control structure
1490		real, intent(in) :: dt !< time step [T ~> s]
1491		type(accel_diag_ptrs), target, intent(inout) :: Accel_diag !< points to momentum equation terms for
1492		!! budget analysis
1493		type(cont_diag_ptrs), target, intent(inout) :: Cont_diag !< points to terms in continuity equation
1494		type(ocean_internal_state), intent(inout) :: MIS !< "MOM6 internal state" used to pass
1495		!! diagnostic pointers
1496		type(VarMix_CS), intent(inout) :: VarMix !< points to spatially variable viscosities
1497		type(MEKE_type), intent(inout) :: MEKE !< MEKE fields
1498		type(thickness_diffuse_CS), intent(inout) :: thickness_diffuse_CSp !< Pointer to the control structure
1499		!! used for the isopycnal height diffusive transport.
1500		type(ocean_OBC_type), pointer :: OBC !< points to OBC related fields
1501		type(update_OBC_CS), pointer :: update_OBC_CSp !< points to OBC update related fields
1502		type(ALE_CS), pointer :: ALE_CSp !< points to ALE control structure
1503		type(set_visc_CS), target, intent(in) :: set_visc !< set_visc control structure
1504		type(vertvisc_type), intent(inout) :: visc !< vertical viscosities, bottom drag, and related
1505		type(directories), intent(in) :: dirs !< contains directory paths
1506		integer, target, intent(inout) :: ntrunc !< A target for the variable that records
1507		!! the number of times the velocity is
1508		!! truncated (this should be 0).
1509		logical, intent(out) :: calc_dtbt !< If true, recalculate the barotropic time step
1510		type(porous_barrier_type), intent(in) :: pbv !< porous barrier fractional cell metrics
1511		integer, intent(out) :: cont_stencil !< The stencil for thickness
1512		!! from the continuity solver.
1513		integer, intent(out) :: dyn_h_stencil !< The stencil for thickness for the
1514		!! the dynamics based on the continuity
1515		!! solver and Coriolis scheme.
1516
1517		! local variables
1518	1	real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_tmp ! A temporary copy of the layer thicknesses [H ~> m or kg m-2]
1519		character(len=40) :: mdl = "MOM_dynamics_split_RK2" ! This module's name.
1520		! This include declares and sets the variable "version".
1521		# include "version_variable.h"
1522		character(len=48) :: thickness_units, flux_units, eta_rest_name
1523		type(group_pass_type) :: pass_av_h_uvh
1524		logical :: debug_truncations
1525		logical :: read_uv, read_h2
1526		logical :: enable_bugs ! If true, the defaults for recently added bug-fix flags are set to
1527		! recreate the bugs, or if false bugs are only used if actively selected.
1528		logical :: visc_rem_bug ! Stores the value of runtime paramter VISC_REM_BUG.
1529		integer :: cor_stencil
1530
1531		integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz
1532		integer :: IsdB, IedB, JsdB, JedB
1533		integer :: nc ! Number of tidal constituents to be harmonically analyzed
1534	1	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
1535	1	isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
1536	1	IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB
1537
1538	1	if (.not.associated(CS)) call MOM_error(FATAL, &
1539	0	"initialize_dyn_split_RK2 called with an unassociated control structure.")
1540	1	if (CS%module_is_initialized) then
1541		call MOM_error(WARNING, "initialize_dyn_split_RK2 called with a control "// &
1542	0	"structure that has already been initialized.")
1543	0	return
1544		endif
1545	1	CS%module_is_initialized = .true.
1546
1547	1	CS%diag => diag
1548
1549	1	call log_version(param_file, mdl, version, "")
1550		call get_param(param_file, mdl, "TIDES", CS%use_tides, &
1551	1	"If true, apply tidal momentum forcing.", default=.false.)
1552		call get_param(param_file, mdl, "CALCULATE_SAL", CS%calculate_SAL, &
1553	1	"If true, calculate self-attraction and loading.", default=CS%use_tides)
1554		call get_param(param_file, mdl, "USE_HA", CS%use_HA, &
1555	1	"If true, perform inline harmonic analysis.", default=.false.)
1556		call get_param(param_file, mdl, "HA_N_CONST", nc, &
1557		"Number of tidal constituents to be harmonically analyzed.", &
1558	1	default=0, do_not_log=.not.CS%use_HA)
1559	1	if (nc<=0) CS%use_HA = .false.
1560		call get_param(param_file, mdl, "BE", CS%be, &
1561		"If SPLIT is true, BE determines the relative weighting "//&
1562		"of a 2nd-order Runga-Kutta baroclinic time stepping "//&
1563		"scheme (0.5) and a backward Euler scheme (1) that is "//&
1564		"used for the Coriolis and inertial terms. BE may be "//&
1565		"from 0.5 to 1, but instability may occur near 0.5. "//&
1566		"BE is also applicable if SPLIT is false and USE_RK2 "//&
1567	1	"is true.", units="nondim", default=0.6)
1568		call get_param(param_file, mdl, "BEGW", CS%begw, &
1569		"If SPLIT is true, BEGW is a number from 0 to 1 that "//&
1570		"controls the extent to which the treatment of gravity "//&
1571		"waves is forward-backward (0) or simulated backward "//&
1572		"Euler (1). 0 is almost always used. "//&
1573		"If SPLIT is false and USE_RK2 is true, BEGW can be "//&
1574		"between 0 and 0.5 to damp gravity waves.", &
1575	1	units="nondim", default=0.0)
1576		call get_param(param_file, mdl, "SET_DTBT_USE_BT_CONT", CS%dtbt_use_bt_cont, &
1577	1	"If true, use BT_CONT to calculate DTBT if possible.", default=.false.)
1578		call get_param(param_file, mdl, "SPLIT_BOTTOM_STRESS", CS%split_bottom_stress, &
1579		"If true, provide the bottom stress calculated by the "//&
1580	1	"vertical viscosity to the barotropic solver.", default=.false.)
1581		call get_param(param_file, mdl, "BT_USE_LAYER_FLUXES", CS%BT_use_layer_fluxes, &
1582		"If true, use the summed layered fluxes plus an "//&
1583		"adjustment due to the change in the barotropic velocity "//&
1584	1	"in the barotropic continuity equation.", default=.true.)
1585		call get_param(param_file, mdl, "BT_ADJ_CORR_MASS_SRC", CS%BT_adj_corr_mass_src, &
1586		"If true, recalculates the barotropic mass source after "//&
1587		"predictor step. This should make little difference in the "//&
1588		"deep ocean but appears to help for vanished layers. If false, "//&
1589	1	"uses the same mass source as from the predictor step.", default=.true.)
1590		call get_param(param_file, mdl, "STORE_CORIOLIS_ACCEL", CS%store_CAu, &
1591		"If true, calculate the Coriolis accelerations at the end of each "//&
1592		"timestep for use in the predictor step of the next split RK2 timestep.", &
1593	1	default=.true.)
1594		call get_param(param_file, mdl, "FPMIX", CS%fpmix, &
1595		"If true, add non-local momentum flux increments and diffuse down the Eulerian gradient.", &
1596	1	default=.false.)
1597	1	if (CS%fpmix) then
1598		call get_param(param_file, "MOM", "CEMP_NL", CS%Cemp_NL, &
1599		"Empirical coefficient of non-local momentum mixing.", &
1600	0	units="nondim", default=3.6)
1601		endif
1602		call get_param(param_file, mdl, "REMAP_AUXILIARY_VARS", CS%remap_aux, &
1603		"If true, apply ALE remapping to all of the auxiliary 3-dimensional "//&
1604		"variables that are needed to reproduce across restarts, similarly to "//&
1605		"what is already being done with the primary state variables. "//&
1606	1	"The default should be changed to true.", default=.false., do_not_log=.true.)
1607	1	if (CS%remap_aux .and. .not.CS%store_CAu) call MOM_error(FATAL, &
1608	0	"REMAP_AUXILIARY_VARS requires that STORE_CORIOLIS_ACCEL = True.")
1609		call get_param(param_file, mdl, "DEBUG", CS%debug, &
1610		"If true, write out verbose debugging data.", &
1611	1	default=.false., debuggingParam=.true.)
1612		call get_param(param_file, mdl, "OBC_DEBUGGING_TESTS", CS%debug_OBC, &
1613		"If true, do additional calls resetting certain values to help verify the "//&
1614		"correctness of the open boundary condition code.", &
1615	1	default=.false., old_name="DEBUG_OBC", debuggingParam=.true., do_not_log=.true.)
1616		call get_param(param_file, mdl, "DEBUG_TRUNCATIONS", debug_truncations, &
1617	1	default=.false.)
1618		call get_param(param_file, mdl, "ENABLE_BUGS_BY_DEFAULT", enable_bugs, &
1619	1	default=.true., do_not_log=.true.) ! This is logged from MOM.F90.
1620		call get_param(param_file, mdl, "VISC_REM_BUG", visc_rem_bug, &
1621		"If true, visc_rem_[uv] in split mode is incorrectly calculated or accounted "//&
1622		"for in two places. This parameter controls the defaults of two individual "//&
1623		"flags, VISC_REM_TIMESTEP_BUG in MOM_dynamics_split_RK2(b) and "//&
1624	1	"VISC_REM_BT_WEIGHT_BUG in MOM_barotropic.", default=.false.)
1625		call get_param(param_file, mdl, "VISC_REM_TIMESTEP_BUG", CS%visc_rem_dt_bug, &
1626		"If true, recover a bug that uses dt_pred rather than dt in "//&
1627		"vertvisc_remnant() at the end of predictor stage for the following "//&
1628		"continuity() and btstep() calls in the corrector step. Default of this flag "//&
1629	1	"is set by VISC_REM_BUG", default=visc_rem_bug)
1630
1631	8841	allocate(CS%taux_bot(IsdB:IedB,jsd:jed), source=0.0)
1632	8902	allocate(CS%tauy_bot(isd:ied,JsdB:JedB), source=0.0)
1633		!$omp target enter data map(to: CS%taux_bot, CS%tauy_bot)
1634
1635	8841	ALLOC_(CS%uhbt(IsdB:IedB,jsd:jed)) ; CS%uhbt(:,:) = 0.0
1636	8902	ALLOC_(CS%vhbt(isd:ied,JsdB:JedB)) ; CS%vhbt(:,:) = 0.0
1637		!$omp target enter data map(alloc: CS%uhbt, CS%vhbt)
1638	663076	ALLOC_(CS%visc_rem_u(IsdB:IedB,jsd:jed,nz)) ; CS%visc_rem_u(:,:,:) = 0.0
1639	667651	ALLOC_(CS%visc_rem_v(isd:ied,JsdB:JedB,nz)) ; CS%visc_rem_v(:,:,:) = 0.0
1640		!$omp target enter data map(alloc: CS%visc_rem_u, CS%visc_rem_v)
1641	8773	ALLOC_(CS%eta_PF(isd:ied,jsd:jed)) ; CS%eta_PF(:,:) = 0.0
1642	657976	ALLOC_(CS%pbce(isd:ied,jsd:jed,nz)) ; CS%pbce(:,:,:) = 0.0
1643		!$omp target enter data map(to: CS%pbce, CS%eta_PF)
1644	663076	ALLOC_(CS%u_accel_bt(IsdB:IedB,jsd:jed,nz)) ; CS%u_accel_bt(:,:,:) = 0.0
1645	667651	ALLOC_(CS%v_accel_bt(isd:ied,JsdB:JedB,nz)) ; CS%v_accel_bt(:,:,:) = 0.0
1646		!$omp target enter data map(alloc: CS%u_accel_bt, CS%v_accel_bt)
1647	663076	ALLOC_(CS%PFu_Stokes(IsdB:IedB,jsd:jed,nz)) ; CS%PFu_Stokes(:,:,:) = 0.0
1648	667651	ALLOC_(CS%PFv_Stokes(isd:ied,JsdB:JedB,nz)) ; CS%PFv_Stokes(:,:,:) = 0.0
1649
1650	1	MIS%diffu => CS%diffu
1651	1	MIS%diffv => CS%diffv
1652	1	MIS%PFu => CS%PFu
1653	1	MIS%PFv => CS%PFv
1654	1	MIS%CAu => CS%CAu
1655	1	MIS%CAv => CS%CAv
1656	1	MIS%pbce => CS%pbce
1657	1	MIS%u_accel_bt => CS%u_accel_bt
1658	1	MIS%v_accel_bt => CS%v_accel_bt
1659	1	MIS%u_av => CS%u_av
1660	1	MIS%v_av => CS%v_av
1661
1662	1	CS%ADp => Accel_diag
1663	1	CS%CDp => Cont_diag
1664	1	Accel_diag%diffu => CS%diffu
1665	1	Accel_diag%diffv => CS%diffv
1666	1	Accel_diag%PFu => CS%PFu
1667	1	Accel_diag%PFv => CS%PFv
1668	1	Accel_diag%CAu => CS%CAu
1669	1	Accel_diag%CAv => CS%CAv
1670	1	Accel_diag%u_accel_bt => CS%u_accel_bt
1671	1	Accel_diag%v_accel_bt => CS%v_accel_bt
1672
1673	1	allocate(CS%AD_pred)
1674	1	CS%AD_pred%diffu => CS%diffu
1675	1	CS%AD_pred%diffv => CS%diffv
1676	1	CS%AD_pred%PFu => CS%PFu
1677	1	CS%AD_pred%PFv => CS%PFv
1678	1	CS%AD_pred%CAu => CS%CAu_pred
1679	1	CS%AD_pred%CAv => CS%CAv_pred
1680	1	CS%AD_pred%u_accel_bt => CS%u_accel_bt
1681	1	CS%AD_pred%v_accel_bt => CS%v_accel_bt
1682
1683		! Accel_diag%pbce => CS%pbce
1684		! Accel_diag%u_accel_bt => CS%u_accel_bt ; Accel_diag%v_accel_bt => CS%v_accel_bt
1685		! Accel_diag%u_av => CS%u_av ; Accel_diag%v_av => CS%v_av
1686
1687	1	id_clock_pass_init = cpu_clock_id('(Ocean init message passing)', grain=CLOCK_ROUTINE)
1688
1689		!$omp target enter data map(alloc: CS%continuity_CSp)
1690	1	call continuity_init(Time, G, GV, US, param_file, diag, CS%continuity_CSp, CS%OBC)
1691	1	cont_stencil = continuity_stencil(CS%continuity_CSp)
1692	1	call CoriolisAdv_init(Time, G, GV, US, param_file, diag, CS%ADp, CS%CoriolisAdv)
1693	1	cor_stencil = CoriolisAdv_stencil(CS%CoriolisAdv)
1694	1	dyn_h_stencil = max(cont_stencil, CoriolisAdv_stencil(CS%CoriolisAdv))
1695	1	if (CS%calculate_SAL) call SAL_init(h, tv, G, GV, US, param_file, CS%SAL_CSp, restart_CS)
1696	1	if (CS%use_tides) call tidal_forcing_init(Time, G, US, param_file, CS%tides_CSp)
1697	1	if (CS%use_HA) then
1698	0	call HA_init(Time, US, param_file, nc, CS%HA_CSp)
1699	0	HA_CSp => CS%HA_CSp
1700		else
1701	1	HA_CSp => NULL()
1702		endif
1703
1704		!$omp target enter data map(alloc: CS%PressureForce_CSp)
1705		call PressureForce_init(Time, G, GV, US, param_file, diag, CS%PressureForce_CSp, CS%ADp, &
1706	1	CS%SAL_CSp, CS%tides_CSp)
1707
1708		!$omp target enter data map(alloc: CS%hor_visc)
1709	1	call hor_visc_init(Time, G, GV, US, param_file, diag, CS%hor_visc, ADp=CS%ADp)
1710
1711	1	allocate(CS%vertvisc_CSp)
1712		!$omp target enter data map(alloc: CS%vertvisc_CSp)
1713		call vertvisc_init(MIS, Time, G, GV, US, param_file, diag, CS%ADp, dirs, &
1714	1	ntrunc, CS%vertvisc_CSp, CS%fpmix)
1715	1	CS%set_visc_CSp => set_visc
1716	1	call updateCFLtruncationValue(Time, CS%vertvisc_CSp, US, activate=is_new_run(restart_CS) )
1717
1718	1	if (associated(ALE_CSp)) CS%ALE_CSp => ALE_CSp
1719	1	if (associated(OBC)) then
1720	0	CS%OBC => OBC
1721	0	if (OBC%ramp) call update_OBC_ramp(Time, CS%OBC, US, activate=is_new_run(restart_CS) )
1722		endif
1723	1	if (associated(update_OBC_CSp)) CS%update_OBC_CSp => update_OBC_CSp
1724
1725	1	eta_rest_name = "sfc" ; if (.not.GV%Boussinesq) eta_rest_name = "p_bot"
1726	1	if (.not. query_initialized(CS%eta, trim(eta_rest_name), restart_CS)) then
1727		! Estimate eta based on the layer thicknesses - h. With the Boussinesq
1728		! approximation, eta is the free surface height anomaly, while without it
1729		! eta is the mass of ocean per unit area. eta always has the same
1730		! dimensions as h, either m or kg m-3.
1731		! CS%eta(:,:) = 0.0 already from initialization.
1732	1	if (GV%Boussinesq) then
1733	7261	do j=js,je ; do i=is,ie ; CS%eta(i,j) = -GV%Z_to_H * G%bathyT(i,j) ; enddo ; enddo
1734		endif
1735	544576	do k=1,nz ; do j=js,je ; do i=is,ie
1736	544500	CS%eta(i,j) = CS%eta(i,j) + h(i,j,k)
1737		enddo ; enddo ; enddo
1738	1	call set_initialized(CS%eta, trim(eta_rest_name), restart_CS)
1739		!$omp target update to(CS%eta)
1740		endif
1741		! Copy eta into an output array.
1742	7261	do j=js,je ; do i=is,ie ; eta(i,j) = CS%eta(i,j) ; enddo ; enddo
1743
1744		!$omp target enter data map (alloc: CS%barotropic_CSp)
1745		call barotropic_init(u, v, h, Time, G, GV, US, param_file, diag, &
1746		CS%barotropic_CSp, restart_CS, calc_dtbt, CS%BT_cont, &
1747	1	CS%OBC, CS%SAL_CSp, HA_CSp)
1748
1749	1	if (.not. query_initialized(CS%diffu, "diffu", restart_CS) .or. &
1750		.not. query_initialized(CS%diffv, "diffv", restart_CS)) then
1751		!$omp target update to(u, v, h, uh, vh)
1752		call horizontal_viscosity(u, v, h, uh, vh, CS%diffu, CS%diffv, MEKE, VarMix, G, GV, US, CS%hor_visc, &
1753		hor_visc_nkblock(CS%hor_visc, GV), tv, dt, OBC=CS%OBC, BT=CS%barotropic_CSp, &
1754	1	TD=thickness_diffuse_CSp, hu_cont=CS%BT_cont%h_u, hv_cont=CS%BT_cont%h_v)
1755		!$omp target update from(CS%diffu, CS%diffv)
1756	1	call set_initialized(CS%diffu, "diffu", restart_CS)
1757	1	call set_initialized(CS%diffv, "diffv", restart_CS)
1758		endif
1759
1760	1	if (.not. query_initialized(CS%u_av, "u2", restart_CS) .or. &
1761		.not. query_initialized(CS%v_av, "v2", restart_CS)) then
1762	663076	do k=1,nz ; do j=jsd,jed ; do I=IsdB,IedB ; CS%u_av(I,j,k) = u(I,j,k) ; enddo ; enddo ; enddo
1763	667651	do k=1,nz ; do J=JsdB,JedB ; do i=isd,ied ; CS%v_av(i,J,k) = v(i,J,k) ; enddo ; enddo ; enddo
1764	1	call set_initialized(CS%u_av, "u2", restart_CS)
1765	1	call set_initialized(CS%v_av, "v2", restart_CS)
1766		endif
1767
1768		!$omp target enter data map(alloc: h_tmp )
1769	1	if (CS%store_CAu) then
1770	1	if (query_initialized(CS%CAu_pred, "CAu", restart_CS) .and. &
1771		query_initialized(CS%CAv_pred, "CAv", restart_CS)) then
1772	0	CS%CAu_pred_stored = .true.
1773		else
1774		call only_read_from_restarts(uh, vh, 'uh', 'vh', G, restart_CS, stagger=CGRID_NE, &
1775		filename=dirs%input_filename, directory=dirs%restart_input_dir, &
1776	1	success=read_uv, scale=US%m_to_L*2US%T_to_s/GV%H_to_mks)
1777		call only_read_from_restarts('h2', CS%h_av, G, restart_CS, &
1778		filename=dirs%input_filename, directory=dirs%restart_input_dir, &
1779	1	success=read_h2, scale=1.0/GV%H_to_mks)
1780	1	if (read_uv .and. read_h2) then
1781	0	call pass_var(CS%h_av, G%Domain, clock=id_clock_pass_init)
1782		else
1783	1	do concurrent (k=1:nz, j=jsd:jed, i=isd:ied)
1784	661633	h_tmp(i,j,k) = h(i,j,k)
1785		enddo
1786		!$omp target update to(CS%u_av, CS%v_av, CS%h_av, uh, vh)
1787	1	call continuity(CS%u_av, CS%v_av, h, h_tmp, uh, vh, dt, G, GV, US, CS%continuity_CSp, CS%OBC, pbv)
1788		!$omp target update from(CS%u_av, CS%v_av, CS%h_av, uh, vh)
1789
1790		!$omp target update from(h_tmp)
1791	1	call pass_var(h_tmp, G%Domain, clock=id_clock_pass_init)
1792		!$omp target update to(h_tmp)
1793
1794	1	do concurrent (k=1:nz, j=jsd:jed, i=isd:ied)
1795	661633	CS%h_av(i,j,k) = 0.5*(h(i,j,k) + h_tmp(i,j,k))
1796		enddo
1797		endif
1798	1	call pass_vector(CS%u_av, CS%v_av, G%Domain, halo=cor_stencil, clock=id_clock_pass_init, complete=.false.)
1799	1	call pass_vector(uh, vh, G%Domain, halo=cor_stencil, clock=id_clock_pass_init, complete=.true.)
1800		!$omp target update to(CS%u_av, CS%v_av, CS%h_av, uh, vh)
1801		call CorAdCalc(CS%u_av, CS%v_av, CS%h_av, uh, vh, CS%CAu_pred, CS%CAv_pred, CS%OBC, CS%ADp, &
1802	1	G, GV, US, CS%CoriolisAdv, pbv) !, Waves=Waves)
1803		!$omp target update from(CS%CAu_pred, CS%CAv_pred)
1804	1	CS%CAu_pred_stored = .true.
1805		endif
1806		else
1807	0	CS%CAu_pred_stored = .false.
1808		! This call is just here to initialize uh and vh.
1809	0	if (.not. query_initialized(uh, "uh", restart_CS) .or. &
1810		.not. query_initialized(vh, "vh", restart_CS)) then
1811	0	do k=1,nz ; do j=jsd,jed ; do i=isd,ied ; h_tmp(i,j,k) = h(i,j,k) ; enddo ; enddo ; enddo
1812	0	call continuity(u, v, h, h_tmp, uh, vh, dt, G, GV, US, CS%continuity_CSp, CS%OBC, pbv)
1813	0	call pass_var(h_tmp, G%Domain, clock=id_clock_pass_init)
1814	0	do k=1,nz ; do j=jsd,jed ; do i=isd,ied
1815	0	CS%h_av(i,j,k) = 0.5*(h(i,j,k) + h_tmp(i,j,k))
1816		enddo ; enddo ; enddo
1817	0	call set_initialized(uh, "uh", restart_CS)
1818	0	call set_initialized(vh, "vh", restart_CS)
1819	0	call set_initialized(CS%h_av, "h2", restart_CS)
1820		! Try reading the CAu and CAv fields from the restart file, in case this restart file is
1821		! using a newer format.
1822		call only_read_from_restarts(CS%CAu_pred, CS%CAv_pred, "CAu", "CAv", G, restart_CS, &
1823		stagger=CGRID_NE, filename=dirs%input_filename, directory=dirs%restart_input_dir, &
1824	0	success=read_uv, scale=US%m_s_to_L_T*US%T_to_s)
1825	0	CS%CAu_pred_stored = read_uv
1826		else
1827	0	if (.not. query_initialized(CS%h_av, "h2", restart_CS)) then
1828	0	CS%h_av(:,:,:) = h(:,:,:)
1829	0	call set_initialized(CS%h_av, "h2", restart_CS)
1830		endif
1831		endif
1832		endif
1833		!$omp target exit data map(delete: h_tmp )
1834	1	call cpu_clock_begin(id_clock_pass_init)
1835	1	call create_group_pass(pass_av_h_uvh, CS%u_av, CS%v_av, G%Domain, halo=2)
1836	1	if (CS%CAu_pred_stored) then
1837	1	call create_group_pass(pass_av_h_uvh, CS%CAu_pred, CS%CAv_pred, G%Domain, halo=2)
1838		else
1839	0	call create_group_pass(pass_av_h_uvh, CS%h_av, G%Domain, halo=2)
1840	0	call create_group_pass(pass_av_h_uvh, uh, vh, G%Domain, halo=2)
1841		endif
1842	1	call do_group_pass(pass_av_h_uvh, G%Domain)
1843	1	call cpu_clock_end(id_clock_pass_init)
1844
1845	1	flux_units = get_flux_units(GV)
1846	1	thickness_units = get_thickness_units(GV)
1847		CS%id_uh = register_diag_field('ocean_model', 'uh', diag%axesCuL, Time, &
1848		'Zonal Thickness Flux', flux_units, conversion=GV%H_to_MKSUS%L_to_m2US%s_to_T, &
1849	1	y_cell_method='sum', v_extensive=.true.)
1850		CS%id_vh = register_diag_field('ocean_model', 'vh', diag%axesCvL, Time, &
1851		'Meridional Thickness Flux', flux_units, conversion=GV%H_to_MKSUS%L_to_m2US%s_to_T, &
1852	1	x_cell_method='sum', v_extensive=.true.)
1853
1854		CS%id_CAu = register_diag_field('ocean_model', 'CAu', diag%axesCuL, Time, &
1855	1	'Zonal Coriolis and Advective Acceleration', 'm s-2', conversion=US%L_T2_to_m_s2)
1856		CS%id_CAv = register_diag_field('ocean_model', 'CAv', diag%axesCvL, Time, &
1857	1	'Meridional Coriolis and Advective Acceleration', 'm s-2', conversion=US%L_T2_to_m_s2)
1858		CS%id_PFu = register_diag_field('ocean_model', 'PFu', diag%axesCuL, Time, &
1859	1	'Zonal Pressure Force Acceleration', 'm s-2', conversion=US%L_T2_to_m_s2)
1860		CS%id_PFv = register_diag_field('ocean_model', 'PFv', diag%axesCvL, Time, &
1861	1	'Meridional Pressure Force Acceleration', 'm s-2', conversion=US%L_T2_to_m_s2)
1862		CS%id_ueffA = register_diag_field('ocean_model', 'ueffA', diag%axesCuL, Time, &
1863		'Effective U-Face Area', 'm^2', conversion=GV%H_to_m*US%L_to_m, &
1864	1	y_cell_method='sum', v_extensive=.true.)
1865		CS%id_veffA = register_diag_field('ocean_model', 'veffA', diag%axesCvL, Time, &
1866		'Effective V-Face Area', 'm^2', conversion=GV%H_to_m*US%L_to_m, &
1867	1	x_cell_method='sum', v_extensive=.true.)
1868	1	if (GV%Boussinesq) then
1869		CS%id_deta_dt = register_diag_field('ocean_model', 'deta_dt', diag%axesT1, Time, &
1870	1	'Barotropic SSH tendency due to dynamics', trim(thickness_units)//' s-1', conversion=GV%H_to_MKS*US%s_to_T)
1871		else
1872		CS%id_deta_dt = register_diag_field('ocean_model', 'deta_dt', diag%axesT1, Time, &
1873		'Barotropic column-mass tendency due to dynamics', trim(thickness_units)//' s-1', &
1874	0	conversion=GV%H_to_mks*US%s_to_T)
1875		endif
1876
1877		!CS%id_hf_PFu = register_diag_field('ocean_model', 'hf_PFu', diag%axesCuL, Time, &
1878		! 'Fractional Thickness-weighted Zonal Pressure Force Acceleration', &
1879		! 'm s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)
1880		!if (CS%id_hf_PFu > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_u,IsdB,IedB,jsd,jed,nz)
1881
1882		!CS%id_hf_PFv = register_diag_field('ocean_model', 'hf_PFv', diag%axesCvL, Time, &
1883		! 'Fractional Thickness-weighted Meridional Pressure Force Acceleration', &
1884		! 'm s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)
1885		!if (CS%id_hf_PFv > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz)
1886
1887		!CS%id_hf_CAu = register_diag_field('ocean_model', 'hf_CAu', diag%axesCuL, Time, &
1888		! 'Fractional Thickness-weighted Zonal Coriolis and Advective Acceleration', &
1889		! 'm s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)
1890		!if (CS%id_hf_CAu > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_u,IsdB,IedB,jsd,jed,nz)
1891
1892		!CS%id_hf_CAv = register_diag_field('ocean_model', 'hf_CAv', diag%axesCvL, Time, &
1893		! 'Fractional Thickness-weighted Meridional Coriolis and Advective Acceleration', &
1894		! 'm s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)
1895		!if (CS%id_hf_CAv > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz)
1896
1897		CS%id_hf_PFu_2d = register_diag_field('ocean_model', 'hf_PFu_2d', diag%axesCu1, Time, &
1898		'Depth-sum Fractional Thickness-weighted Zonal Pressure Force Acceleration', &
1899	1	'm s-2', conversion=US%L_T2_to_m_s2)
1900	1	if (CS%id_hf_PFu_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_u,IsdB,IedB,jsd,jed,nz)
1901
1902		CS%id_hf_PFv_2d = register_diag_field('ocean_model', 'hf_PFv_2d', diag%axesCv1, Time, &
1903		'Depth-sum Fractional Thickness-weighted Meridional Pressure Force Acceleration', &
1904	1	'm s-2', conversion=US%L_T2_to_m_s2)
1905	1	if (CS%id_hf_PFv_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz)
1906
1907		CS%id_h_PFu = register_diag_field('ocean_model', 'h_PFu', diag%axesCuL, Time, &
1908		'Thickness Multiplied Zonal Pressure Force Acceleration', &
1909	1	'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
1910	1	if (CS%id_h_PFu > 0) call safe_alloc_ptr(CS%ADp%diag_hu,IsdB,IedB,jsd,jed,nz)
1911
1912		CS%id_h_PFv = register_diag_field('ocean_model', 'h_PFv', diag%axesCvL, Time, &
1913		'Thickness Multiplied Meridional Pressure Force Acceleration', &
1914	1	'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
1915	1	if (CS%id_h_PFv > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,JsdB,JedB,nz)
1916
1917		CS%id_intz_PFu_2d = register_diag_field('ocean_model', 'intz_PFu_2d', diag%axesCu1, Time, &
1918		'Depth-integral of Zonal Pressure Force Acceleration', &
1919	1	'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
1920	1	if (CS%id_intz_PFu_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hu,IsdB,IedB,jsd,jed,nz)
1921
1922		CS%id_intz_PFv_2d = register_diag_field('ocean_model', 'intz_PFv_2d', diag%axesCv1, Time, &
1923		'Depth-integral of Meridional Pressure Force Acceleration', &
1924	1	'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
1925	1	if (CS%id_intz_PFv_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,JsdB,JedB,nz)
1926
1927		CS%id_hf_CAu_2d = register_diag_field('ocean_model', 'hf_CAu_2d', diag%axesCu1, Time, &
1928		'Depth-sum Fractional Thickness-weighted Zonal Coriolis and Advective Acceleration', &
1929	1	'm s-2', conversion=US%L_T2_to_m_s2)
1930	1	if (CS%id_hf_CAu_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_u,IsdB,IedB,jsd,jed,nz)
1931
1932		CS%id_hf_CAv_2d = register_diag_field('ocean_model', 'hf_CAv_2d', diag%axesCv1, Time, &
1933		'Depth-sum Fractional Thickness-weighted Meridional Coriolis and Advective Acceleration', &
1934	1	'm s-2', conversion=US%L_T2_to_m_s2)
1935	1	if (CS%id_hf_CAv_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz)
1936
1937		CS%id_h_CAu = register_diag_field('ocean_model', 'h_CAu', diag%axesCuL, Time, &
1938		'Thickness Multiplied Zonal Coriolis and Advective Acceleration', &
1939	1	'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
1940	1	if (CS%id_h_CAu > 0) call safe_alloc_ptr(CS%ADp%diag_hu,IsdB,IedB,jsd,jed,nz)
1941
1942		CS%id_h_CAv = register_diag_field('ocean_model', 'h_CAv', diag%axesCvL, Time, &
1943		'Thickness Multiplied Meridional Coriolis and Advective Acceleration', &
1944	1	'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
1945	1	if (CS%id_h_CAv > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,JsdB,JedB,nz)
1946
1947		CS%id_intz_CAu_2d = register_diag_field('ocean_model', 'intz_CAu_2d', diag%axesCu1, Time, &
1948		'Depth-integral of Zonal Coriolis and Advective Acceleration', &
1949	1	'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
1950	1	if (CS%id_intz_CAu_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hu,IsdB,IedB,jsd,jed,nz)
1951
1952		CS%id_intz_CAv_2d = register_diag_field('ocean_model', 'intz_CAv_2d', diag%axesCv1, Time, &
1953		'Depth-integral of Meridional Coriolis and Advective Acceleration', &
1954	1	'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
1955	1	if (CS%id_intz_CAv_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,JsdB,JedB,nz)
1956
1957		CS%id_uav = register_diag_field('ocean_model', 'uav', diag%axesCuL, Time, &
1958	1	'Barotropic-step Averaged Zonal Velocity', 'm s-1', conversion=US%L_T_to_m_s)
1959		CS%id_vav = register_diag_field('ocean_model', 'vav', diag%axesCvL, Time, &
1960	1	'Barotropic-step Averaged Meridional Velocity', 'm s-1', conversion=US%L_T_to_m_s)
1961
1962		CS%id_u_BT_accel = register_diag_field('ocean_model', 'u_BT_accel', diag%axesCuL, Time, &
1963	1	'Barotropic Anomaly Zonal Acceleration', 'm s-2', conversion=US%L_T2_to_m_s2)
1964		CS%id_v_BT_accel = register_diag_field('ocean_model', 'v_BT_accel', diag%axesCvL, Time, &
1965	1	'Barotropic Anomaly Meridional Acceleration', 'm s-2', conversion=US%L_T2_to_m_s2)
1966
1967		!CS%id_hf_u_BT_accel = register_diag_field('ocean_model', 'hf_u_BT_accel', diag%axesCuL, Time, &
1968		! 'Fractional Thickness-weighted Barotropic Anomaly Zonal Acceleration', &
1969		! 'm s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)
1970		!if (CS%id_hf_u_BT_accel > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_u,IsdB,IedB,jsd,jed,nz)
1971
1972		!CS%id_hf_v_BT_accel = register_diag_field('ocean_model', 'hf_v_BT_accel', diag%axesCvL, Time, &
1973		! 'Fractional Thickness-weighted Barotropic Anomaly Meridional Acceleration', &
1974		! 'm s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)
1975		!if (CS%id_hf_v_BT_accel > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz)
1976
1977		CS%id_hf_u_BT_accel_2d = register_diag_field('ocean_model', 'hf_u_BT_accel_2d', diag%axesCu1, Time, &
1978		'Depth-sum Fractional Thickness-weighted Barotropic Anomaly Zonal Acceleration', &
1979	1	'm s-2', conversion=US%L_T2_to_m_s2)
1980	1	if (CS%id_hf_u_BT_accel_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_u,IsdB,IedB,jsd,jed,nz)
1981
1982		CS%id_hf_v_BT_accel_2d = register_diag_field('ocean_model', 'hf_v_BT_accel_2d', diag%axesCv1, Time, &
1983		'Depth-sum Fractional Thickness-weighted Barotropic Anomaly Meridional Acceleration', &
1984	1	'm s-2', conversion=US%L_T2_to_m_s2)
1985	1	if (CS%id_hf_v_BT_accel_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz)
1986
1987		CS%id_h_u_BT_accel = register_diag_field('ocean_model', 'h_u_BT_accel', diag%axesCuL, Time, &
1988		'Thickness Multiplied Barotropic Anomaly Zonal Acceleration', &
1989	1	'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
1990	1	if (CS%id_h_u_BT_accel > 0) call safe_alloc_ptr(CS%ADp%diag_hu,IsdB,IedB,jsd,jed,nz)
1991
1992		CS%id_h_v_BT_accel = register_diag_field('ocean_model', 'h_v_BT_accel', diag%axesCvL, Time, &
1993		'Thickness Multiplied Barotropic Anomaly Meridional Acceleration', &
1994	1	'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
1995	1	if (CS%id_h_v_BT_accel > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,JsdB,JedB,nz)
1996
1997		CS%id_intz_u_BT_accel_2d = register_diag_field('ocean_model', 'intz_u_BT_accel_2d', diag%axesCu1, Time, &
1998		'Depth-integral of Barotropic Anomaly Zonal Acceleration', &
1999	1	'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
2000	1	if (CS%id_intz_u_BT_accel_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hu,IsdB,IedB,jsd,jed,nz)
2001
2002		CS%id_intz_v_BT_accel_2d = register_diag_field('ocean_model', 'intz_v_BT_accel_2d', diag%axesCv1, Time, &
2003		'Depth-integral of Barotropic Anomaly Meridional Acceleration', &
2004	1	'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
2005	1	if (CS%id_intz_v_BT_accel_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,JsdB,JedB,nz)
2006
2007		CS%id_PFu_visc_rem = register_diag_field('ocean_model', 'PFu_visc_rem', diag%axesCuL, Time, &
2008		'Zonal Pressure Force Acceleration multiplied by the viscous remnant', &
2009	1	'm s-2', conversion=US%L_T2_to_m_s2)
2010	1	if (CS%id_PFu_visc_rem > 0) call safe_alloc_ptr(CS%ADp%visc_rem_u,IsdB,IedB,jsd,jed,nz)
2011		CS%id_PFv_visc_rem = register_diag_field('ocean_model', 'PFv_visc_rem', diag%axesCvL, Time, &
2012		'Meridional Pressure Force Acceleration multiplied by the viscous remnant', &
2013	1	'm s-2', conversion=US%L_T2_to_m_s2)
2014	1	if (CS%id_PFv_visc_rem > 0) call safe_alloc_ptr(CS%ADp%visc_rem_v,isd,ied,JsdB,JedB,nz)
2015
2016		CS%id_CAu_visc_rem = register_diag_field('ocean_model', 'CAu_visc_rem', diag%axesCuL, Time, &
2017		'Zonal Coriolis and Advective Acceleration multiplied by the viscous remnant', &
2018	1	'm s-2', conversion=US%L_T2_to_m_s2)
2019	1	if (CS%id_CAu_visc_rem > 0) call safe_alloc_ptr(CS%ADp%visc_rem_u,IsdB,IedB,jsd,jed,nz)
2020		CS%id_CAv_visc_rem = register_diag_field('ocean_model', 'CAv_visc_rem', diag%axesCvL, Time, &
2021		'Meridional Coriolis and Advective Acceleration multiplied by the viscous remnant', &
2022	1	'm s-2', conversion=US%L_T2_to_m_s2)
2023	1	if (CS%id_CAv_visc_rem > 0) call safe_alloc_ptr(CS%ADp%visc_rem_v,isd,ied,JsdB,JedB,nz)
2024
2025		CS%id_u_BT_accel_visc_rem = register_diag_field('ocean_model', 'u_BT_accel_visc_rem', diag%axesCuL, Time, &
2026		'Barotropic Anomaly Zonal Acceleration multiplied by the viscous remnant', &
2027	1	'm s-2', conversion=US%L_T2_to_m_s2)
2028	1	if (CS%id_u_BT_accel_visc_rem > 0) call safe_alloc_ptr(CS%ADp%visc_rem_u,IsdB,IedB,jsd,jed,nz)
2029		CS%id_v_BT_accel_visc_rem = register_diag_field('ocean_model', 'v_BT_accel_visc_rem', diag%axesCvL, Time, &
2030		'Barotropic Anomaly Meridional Acceleration multiplied by the viscous remnant', &
2031	1	'm s-2', conversion=US%L_T2_to_m_s2)
2032	1	if (CS%id_v_BT_accel_visc_rem > 0) call safe_alloc_ptr(CS%ADp%visc_rem_v,isd,ied,JsdB,JedB,nz)
2033
2034	1	id_clock_Cor = cpu_clock_id('(Ocean Coriolis & mom advection)', grain=CLOCK_MODULE)
2035	1	id_clock_continuity = cpu_clock_id('(Ocean continuity equation)', grain=CLOCK_MODULE)
2036	1	id_clock_pres = cpu_clock_id('(Ocean pressure force)', grain=CLOCK_MODULE)
2037	1	id_clock_vertvisc = cpu_clock_id('(Ocean vertical viscosity)', grain=CLOCK_MODULE)
2038	1	id_clock_horvisc = cpu_clock_id('(Ocean horizontal viscosity)', grain=CLOCK_MODULE)
2039	1	id_clock_mom_update = cpu_clock_id('(Ocean momentum increments)', grain=CLOCK_MODULE)
2040	1	id_clock_pass = cpu_clock_id('(Ocean message passing)', grain=CLOCK_MODULE)
2041	1	id_clock_btcalc = cpu_clock_id('(Ocean barotropic mode calc)', grain=CLOCK_MODULE)
2042	1	id_clock_btstep = cpu_clock_id('(Ocean barotropic mode stepping)', grain=CLOCK_MODULE)
2043	1	id_clock_btforce = cpu_clock_id('(Ocean barotropic forcing calc)', grain=CLOCK_MODULE)
2044
2045	2	end subroutine initialize_dyn_split_RK2
2046
2047
2048		!> Close the dyn_split_RK2 module
2049	1	subroutine end_dyn_split_RK2(CS)
2050		type(MOM_dyn_split_RK2_CS), pointer :: CS !< module control structure
2051
2052		!$omp target exit data map(delete: CS%barotropic_CSp)
2053	1	call barotropic_end(CS%barotropic_CSp)
2054
2055	1	call vertvisc_end(CS%vertvisc_CSp)
2056	1	deallocate(CS%vertvisc_CSp)
2057
2058	1	call hor_visc_end(CS%hor_visc)
2059		!$omp target exit data map(delete: CS%hor_visc)
2060
2061	1	if (CS%calculate_SAL) call SAL_end(CS%SAL_CSp)
2062	1	if (CS%use_tides) call tidal_forcing_end(CS%tides_CSp)
2063	1	call CoriolisAdv_end(CS%CoriolisAdv)
2064
2065	1	DEALLOC_(CS%diffu) ; DEALLOC_(CS%diffv)
2066		!$omp target exit data map(delete: CS%diffu, CS%diffv)
2067	1	DEALLOC_(CS%CAu) ; DEALLOC_(CS%CAv)
2068		!$omp target exit data map(delete: CS%CAu, CS%CAv)
2069	1	DEALLOC_(CS%CAu_pred) ; DEALLOC_(CS%CAv_pred)
2070		!$omp target exit data map(delete: CS%CAu_pred, CS%CAv_pred)
2071	1	DEALLOC_(CS%PFu) ; DEALLOC_(CS%PFv)
2072		!$omp target exit data map(delete: CS%PFu, CS%PFv)
2073
2074	1	if (associated(CS%taux_bot)) deallocate(CS%taux_bot)
2075	1	if (associated(CS%tauy_bot)) deallocate(CS%tauy_bot)
2076	1	DEALLOC_(CS%uhbt) ; DEALLOC_(CS%vhbt)
2077	1	DEALLOC_(CS%u_accel_bt) ; DEALLOC_(CS%v_accel_bt)
2078	1	DEALLOC_(CS%visc_rem_u) ; DEALLOC_(CS%visc_rem_v)
2079
2080	1	DEALLOC_(CS%eta) ; DEALLOC_(CS%eta_PF) ; DEALLOC_(CS%pbce)
2081		!$omp target exit data map(delete: CS%eta, CS%eta_PF, CS%pbce)
2082	1	DEALLOC_(CS%h_av) ; DEALLOC_(CS%u_av) ; DEALLOC_(CS%v_av)
2083		!$omp target exit data map(delete: CS%u_av, CS%v_av, CS%h_av)
2084
2085	1	call dealloc_BT_cont_type(CS%BT_cont)
2086	1	deallocate(CS%AD_pred)
2087
2088	1	deallocate(CS)
2089	1	end subroutine end_dyn_split_RK2
2090
2091
2092		!> \namespace mom_dynamics_split_rk2
2093		!!
2094		!! This file time steps the adiabatic dynamic core by splitting
2095		!! between baroclinic and barotropic modes. It uses a pseudo-second order
2096		!! Runge-Kutta time stepping scheme for the baroclinic momentum
2097		!! equation and a forward-backward coupling between the baroclinic
2098		!! momentum and continuity equations. This split time-stepping
2099		!! scheme is described in detail in Hallberg (JCP, 1997). Additional
2100		!! issues related to exact tracer conservation and how to
2101		!! ensure consistency between the barotropic and layered estimates
2102		!! of the free surface height are described in Hallberg and
2103		!! Adcroft (Ocean Modelling, 2009). This was the time stepping code
2104		!! that is used for most GOLD applications, including GFDL's ESM2G
2105		!! Earth system model, and all of the examples provided with the
2106		!! MOM code (although several of these solutions are routinely
2107		!! verified by comparison with the slower unsplit schemes).
2108		!!
2109		!! The subroutine step_MOM_dyn_split_RK2 actually does the time
2110		!! stepping, while register_restarts_dyn_split_RK2 sets the fields
2111		!! that are found in a full restart file with this scheme, and
2112		!! initialize_dyn_split_RK2 initializes the cpu clocks that are
2113		!! used in this module. For largely historical reasons, this module
2114		!! does not have its own control structure, but shares the same
2115		!! control structure with MOM.F90 and the other MOM_dynamics_...
2116		!! modules.
2117
2118	0	end module MOM_dynamics_split_RK2