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	4	subroutine step_MOM_dyn_split_RK2(u_inst, v_inst, h, tv, visc, Time_local, dt, forces, &
303	4	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	8	real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: up ! Predicted zonal velocity [L T-1 ~> m s-1].
350	8	real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: vp ! Predicted meridional velocity [L T-1 ~> m s-1].
351	8	real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: hp ! Predicted thickness [H ~> m or kg m-2].
352	8	real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: dz ! Distance between the interfaces around a layer [Z ~> m]
353	8	real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: ueffA ! Effective Area of U-Faces [H L ~> m2]
354	8	real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: veffA ! Effective Area of V-Faces [H L ~> m2]
355	8	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	8	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	8	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	8	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	8	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	8	real, dimension(SZI_(G),SZJ_(G)) :: SpV_avg ! The column averaged specific volume [R-1 ~> m3 kg-1]
370	8	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	8	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	8	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	8	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	8	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	4	u_av, & ! The zonal velocity time-averaged over a time step [L T-1 ~> m s-1].
403	4	v_av, & ! The meridional velocity time-averaged over a time step [L T-1 ~> m s-1].
404	4	h_av ! The layer thickness time-averaged over a time step [H ~> m or kg m-2].
405
406	8	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	8	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	4	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
425	4	Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
426	4	u_av => CS%u_av ; v_av => CS%v_av ; h_av => CS%h_av ; eta => CS%eta
427
428	4	Idt_bc = 1.0 / dt
429
430	4	sym = G%Domain%symmetric ! switch to include symmetric domain in checksums
431
432	4	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	1480	do concurrent (k=1:nz, j=G%jsd:G%jed, I=G%IsdB:G%IedB)
439	6383704	up(I,j,k) = 0.0
440		enddo
441	1476	do concurrent (k=1:nz, J=G%JsdB:G%JedB, i=G%isd:G%ied)
442	6400260	vp(i,J,k) = 0.0
443		enddo
444	1476	do concurrent (k=1:nz, j=G%jsd:G%jed, i=G%isd:G%ied)
445	6366404	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	4	call updateCFLtruncationValue(Time_local, CS%vertvisc_CSp, US)
451
452	4	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	4	dyn_p_surf = associated(p_surf_begin) .and. associated(p_surf_end)
462	4	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	4	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	4	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	4	BT_cont_BT_thick = .false.
489	4	if (associated(CS%BT_cont)) BT_cont_BT_thick = &
490	4	(allocated(CS%BT_cont%h_u) .and. allocated(CS%BT_cont%h_v))
491
492	4	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	4	cont_stencil = continuity_stencil(CS%continuity_CSp)
499	4	obc_stencil = 2
500	4	cor_stencil = CoriolisAdv_stencil(CS%CoriolisAdv)
501	4	if (associated(CS%OBC)) then
502	0	if (CS%OBC%oblique_BCs_exist_globally) obc_stencil = 3
503		endif
504	4	vel_stencil = max(2, obc_stencil, hor_visc_vel_stencil(CS%hor_visc))
505	4	call cpu_clock_begin(id_clock_pass)
506	4	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	4	To_All+SCALAR_PAIR, CGRID_NE, halo=max(1,cont_stencil))
509	4	call create_group_pass(CS%pass_uvp, up, vp, G%Domain, halo=max(1,cont_stencil))
510	4	call create_group_pass(CS%pass_uv, u_inst, v_inst, G%Domain, halo=max(2,cont_stencil))
511
512	4	call create_group_pass(CS%pass_hp_uv, hp, G%Domain, halo=cor_stencil)
513	4	call create_group_pass(CS%pass_hp_uv, u_av, v_av, G%Domain, halo=max(cor_stencil,vel_stencil))
514	4	call create_group_pass(CS%pass_hp_uv, uh(:,:,:), vh(:,:,:), G%Domain, halo=max(cor_stencil,vel_stencil))
515	4	call create_group_pass(CS%pass_h, h, g%domain, halo=max(cor_stencil,cont_stencil))
516	4	call create_group_pass(CS%pass_av_uvh, u_av, v_av, g%domain, halo=max(cor_stencil,vel_stencil))
517	4	call create_group_pass(CS%pass_av_uvh, uh(:,:,:), vh(:,:,:), G%Domain, halo=max(cor_stencil,vel_stencil))
518
519	4	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	4	if (CS%begw == 0.0) call enable_averages(dt, Time_local, CS%diag)
525	4	call cpu_clock_begin(id_clock_pres)
526
527		call PressureForce(h, tv, CS%PFu, CS%PFv, G, GV, US, CS%PressureForce_CSp, &
528	4	CS%ALE_CSp, CS%ADp, p_surf, CS%pbce, CS%eta_PF)
529
530	4	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	4	Use_Stokes_PGF = present(Waves)
540	4	if (Use_Stokes_PGF) then
541	4	Use_Stokes_PGF = associated(Waves)
542	4	if (Use_Stokes_PGF) Use_Stokes_PGF = Waves%Stokes_PGF
543	4	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	4	call cpu_clock_end(id_clock_pres)
566	4	call disable_averaging(CS%diag)
567	4	if (showCallTree) call callTree_wayPoint("done with PressureForce (step_MOM_dyn_split_RK2)")
568
569	4	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	4	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	4	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	4	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	4	call cpu_clock_begin(id_clock_btforce)
597
598	4	do concurrent (k=1:nz, j=js:je, I=Isq:Ieq)
599	5979608	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	4	do concurrent (k=1:nz, J=Jsq:Jeq, i=is:ie)
602	5996164	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	4	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	4	call cpu_clock_end(id_clock_btforce)
611
612	4	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	4	call cpu_clock_begin(id_clock_vertvisc)
629
630	4	do concurrent (k=1:nz, j=js:je, I=Isq:Ieq)
631	5979608	up(I,j,k) = G%mask2dCu(I,j) * (u_inst(I,j,k) + dt * u_bc_accel(I,j,k))
632		enddo
633	4	do concurrent (k=1:nz, J=Jsq:Jeq, i=is:ie)
634	5996164	vp(i,J,k) = G%mask2dCv(i,J) * (v_inst(i,J,k) + dt * v_bc_accel(i,J,k))
635		enddo
636
637	4	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	4	call set_viscous_ML(u_inst, v_inst, h, tv, forces, visc, dt, G, GV, US, CS%set_visc_CSp)
641	4	call disable_averaging(CS%diag)
642
643	4	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	4	call thickness_to_dz(h, tv, dz, G, GV, US, halo_size=1, do_offload=.true.)
649
650	4	call vertvisc_coef(up, vp, h, dz, forces, visc, tv, dt, G, GV, US, CS%vertvisc_CSp, CS%OBC, VarMix)
651	4	call vertvisc_remnant(visc, CS%visc_rem_u, CS%visc_rem_v, dt, G, GV, US, CS%vertvisc_CSp)
652
653	4	call cpu_clock_end(id_clock_vertvisc)
654	4	if (showCallTree) call callTree_wayPoint("done with vertvisc_coef (step_MOM_dyn_split_RK2)")
655
656	4	call cpu_clock_begin(id_clock_pass)
657	4	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	4	call do_group_pass(CS%pass_eta, G%Domain, omp_offload=.true.)
664	4	call do_group_pass(CS%pass_visc_rem, G%Domain, omp_offload=.true.)
665		endif
666	4	call cpu_clock_end(id_clock_pass)
667
668	4	call cpu_clock_begin(id_clock_btcalc)
669		! Calculate the relative layer weights for determining barotropic quantities.
670	4	if (.not.BT_cont_BT_thick) then
671	0	call btcalc(h, G, GV, CS%barotropic_CSp, OBC=CS%OBC)
672		endif
673	4	call bt_mass_source(h, eta, .true., G, GV, CS%barotropic_CSp)
674
675	277492	SpV_avg(:,:) = 0.0
676	4	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	4	call cpu_clock_end(id_clock_btcalc)
683
684	4	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	4	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	4	if (associated(CS%BT_cont) .or. CS%BT_use_layer_fluxes) then
694	4	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	4	visc_rem_u=CS%visc_rem_u, visc_rem_v=CS%visc_rem_v, BT_cont=CS%BT_cont)
697
698	4	call cpu_clock_end(id_clock_continuity)
699	4	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	4	OBC=CS%OBC)
703
704		endif
705	4	if (showCallTree) call callTree_wayPoint("done with continuity[BT_cont] (step_MOM_dyn_split_RK2)")
706		endif
707
708	4	if (CS%BT_use_layer_fluxes) then
709	4	uh_ptr => uh_in ; vh_ptr => vh_in ; u_ptr => u_inst ; v_ptr => v_inst
710		endif
711
712	4	call cpu_clock_begin(id_clock_btstep)
713	4	if (calc_dtbt) then
714	4	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	4	call set_dtbt(G, GV, US, CS%barotropic_CSp, CS%pbce, eta=eta)
720		endif
721		endif
722
723	4	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	4	eta_PF_start, taux_bot, tauy_bot, uh_ptr, vh_ptr, u_ptr, v_ptr)
730
731	4	if (showCallTree) call callTree_leave("btstep()")
732	4	call cpu_clock_end(id_clock_btstep)
733
734		! up = u + dt_pred*( u_bc_accel + u_accel_bt )
735	4	dt_pred = dt * CS%be
736	4	call cpu_clock_begin(id_clock_mom_update)
737
738	4	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	5996164	(v_bc_accel(i,J,k) + CS%v_accel_bt(i,J,k)))
741		enddo
742
743	4	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	5979608	(u_bc_accel(I,j,k) + CS%u_accel_bt(I,j,k)))
746		enddo
747
748	4	call cpu_clock_end(id_clock_mom_update)
749
750	4	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	4	call cpu_clock_begin(id_clock_vertvisc)
777	4	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	4	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	4	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	4	CS%OBC, VarMix)
803
804	4	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	4	GV, US, CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot, waves=waves)
822
823		endif
824
825	4	if (showCallTree) call callTree_wayPoint("done with vertvisc (step_MOM_dyn_split_RK2)")
826	4	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	4	if (CS%visc_rem_dt_bug) then
834	4	call vertvisc_remnant(visc, CS%visc_rem_u, CS%visc_rem_v, dt_pred, G, GV, US, CS%vertvisc_CSp)
835		else
836	0	call vertvisc_remnant(visc, CS%visc_rem_u, CS%visc_rem_v, dt, G, GV, US, CS%vertvisc_CSp)
837		endif
838
839	4	call cpu_clock_end(id_clock_vertvisc)
840
841	4	call do_group_pass(CS%pass_visc_rem, G%Domain, clock=id_clock_pass, omp_offload=.true.)
842
843	4	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	4	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	4	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	4	u_cor=u_av, v_cor=v_av, BT_cont=CS%BT_cont)
856	4	call cpu_clock_end(id_clock_continuity)
857
858	4	if (showCallTree) call callTree_wayPoint("done with continuity (step_MOM_dyn_split_RK2)")
859
860	4	call do_group_pass(CS%pass_hp_uv, G%Domain, clock=id_clock_pass, omp_offload=.true.)
861
862	4	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	1460	do concurrent (k=1:nz, j=js-cor_stencil:je+cor_stencil, i=is-cor_stencil:ie+cor_stencil)
882	6163252	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	4	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	4	if (CS%BT_adj_corr_mass_src) then
893	4	call cpu_clock_begin(id_clock_btcalc)
894	4	call bt_mass_source(hp, eta_pred, .false., G, GV, CS%barotropic_CSp)
895	4	call cpu_clock_end(id_clock_btcalc)
896		endif
897
898	4	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	4	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	4	OBC=CS%OBC)
944	4	if (showCallTree) call callTree_wayPoint("done with btcalc[BT_cont_BT_thick] (step_MOM_dyn_split_RK2)")
945		endif
946
947	4	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	4	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	4	ADp=CS%ADp, hu_cont=CS%BT_cont%h_u, hv_cont=CS%BT_cont%h_v, STOCH=STOCH)
966
967	4	call cpu_clock_end(id_clock_horvisc)
968	4	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	4	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	4	G, GV, US, CS%CoriolisAdv, pbv, Waves=Waves)
974	4	call cpu_clock_end(id_clock_Cor)
975
976	4	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	4	call cpu_clock_begin(id_clock_btforce)
982
983	4	do concurrent (k=1:nz, j=js:je, I=Isq:Ieq)
984	5979608	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	4	do concurrent (k=1:nz, J=Jsq:Jeq,i=is:ie)
988	5996164	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	4	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	4	call cpu_clock_end(id_clock_btforce)
997
998	4	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	4	call cpu_clock_begin(id_clock_btstep)
1017	4	if (CS%BT_use_layer_fluxes) then
1018	4	uh_ptr => uh ; vh_ptr => vh ; u_ptr => u_av ; v_ptr => v_av
1019		endif
1020
1021	4	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	4	eta_PF_start, taux_bot, tauy_bot, uh_ptr, vh_ptr, u_ptr, v_ptr, etaav=eta_av)
1027
1028	4	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	4	do concurrent (j=js:je, i=is:ie)
1034	260644	eta(i,j) = eta_pred(i,j)
1035		enddo
1036
1037	4	call cpu_clock_end(id_clock_btstep)
1038	4	if (showCallTree) call callTree_leave("btstep()")
1039
1040	4	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	4	call cpu_clock_begin(id_clock_mom_update)
1047	4	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	5979608	(u_bc_accel(I,j,k) + CS%u_accel_bt(I,j,k)))
1050		enddo
1051	4	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	5996164	(v_bc_accel(i,J,k) + CS%v_accel_bt(i,J,k)))
1054		enddo
1055	4	call cpu_clock_end(id_clock_mom_update)
1056
1057	4	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	4	call cpu_clock_begin(id_clock_vertvisc)
1078
1079	4	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	4	call thickness_to_dz(h, tv, dz, G, GV, US, halo_size=1, do_offload=.true.)
1098
1099	4	call vertvisc_coef(u_inst, v_inst, h, dz, forces, visc, tv, dt, G, GV, US, CS%vertvisc_CSp, CS%OBC, VarMix)
1100
1101	4	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	4	CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot, waves=waves)
1113		endif
1114
1115	4	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	4	call vertvisc_remnant(visc, CS%visc_rem_u, CS%visc_rem_v, dt, G, GV, US, CS%vertvisc_CSp)
1122
1123	4	call cpu_clock_end(id_clock_vertvisc)
1124	4	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	1460	do concurrent (k=1:nz, j=js-cor_stencil:je+cor_stencil, i=is-cor_stencil:ie+cor_stencil)
1128	6163252	h_av(i,j,k) = h(i,j,k)
1129		enddo
1130
1131	4	call do_group_pass(CS%pass_visc_rem, G%Domain, clock=id_clock_pass, omp_offload=.true.)
1132
1133	4	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	4	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	4	call cpu_clock_begin(id_clock_continuity)
1144	1476	do concurrent (k=1:nz, j=G%jsd:G%jed, i=G%isd:G%ied)
1145	6366404	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	4	u_cor=u_av, v_cor=v_av)
1151	4	call cpu_clock_end(id_clock_continuity)
1152
1153	4	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	4	call diag_update_remap_grids(CS%diag)
1159	4	if (showCallTree) call callTree_wayPoint("done with continuity (step_MOM_dyn_split_RK2)")
1160
1161	4	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	4	call do_group_pass(CS%pass_av_uvh, G%domain, clock=id_clock_pass, omp_offload=.true.)
1166		endif
1167
1168	4	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	1460	do concurrent (k=1:nz, j=js-cor_stencil:je+cor_stencil, i=is-cor_stencil:ie+cor_stencil)
1177	6163252	h_av(i,j,k) = 0.5*(h_av(i,j,k) + h(i,j,k))
1178		enddo
1179
1180	4	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	1464	do concurrent (k=1:nz, j=js-cor_stencil:je+cor_stencil, I=Isq-cor_stencil:Ieq+cor_stencil)
1186	6180184	uhtr(I,j,k) = uhtr(I,j,k) + uh(I,j,k) * dt
1187		enddo
1188
1189	1460	do concurrent (k=1:nz, J=Jsq-cor_stencil:Jeq+cor_stencil, i=is-cor_stencil:ie+cor_stencil)
1190	6196740	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	4	if (associated(CS%OBC)) then
1197	0	call update_segment_thickness_reservoirs(G, GV, uhtr, vhtr, h, CS%OBC)
1198		endif
1199
1200	4	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	4	call cpu_clock_begin(id_clock_Cor)
1204	4	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	4	G, GV, US, CS%CoriolisAdv, pbv, Waves=Waves)
1208	4	CS%CAu_pred_stored = .true.
1209	4	call enable_averages(dt, Time_local, CS%diag) ! Reenable the averaging
1210	4	call cpu_clock_end(id_clock_Cor)
1211	4	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	4	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	4	if (CS%id_PFu > 0) call post_data(CS%id_PFu, CS%PFu, CS%diag)
1224	4	if (CS%id_PFv > 0) call post_data(CS%id_PFv, CS%PFv, CS%diag)
1225	4	if (CS%id_CAu > 0) call post_data(CS%id_CAu, CS%CAu, CS%diag)
1226	4	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	4	if (CS%id_uh > 0) call post_data(CS%id_uh, uh, CS%diag)
1230	4	if (CS%id_vh > 0) call post_data(CS%id_vh, vh, CS%diag)
1231	4	if (CS%id_uav > 0) call post_data(CS%id_uav, u_av, CS%diag)
1232	4	if (CS%id_vav > 0) call post_data(CS%id_vav, v_av, CS%diag)
1233	4	if (CS%id_u_BT_accel > 0) call post_data(CS%id_u_BT_accel, CS%u_accel_bt, CS%diag)
1234	4	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	4	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	4	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	4	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	4	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	4	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	4	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	4	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	4	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	4	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	4	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	4	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	4	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	4	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	4	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	4	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	4	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	4	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	4	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	4	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	4	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	4	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	4	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	4	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	4	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	4	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	4	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	4	if (CS%id_deta_dt > 0) call post_data(CS%id_deta_dt, deta_dt, CS%diag)
1313
1314	4	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	4	if (showCallTree) call callTree_leave("step_MOM_dyn_split_RK2()")
1323
1324	4	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	1530343	ALLOC_(CS%diffu(IsdB:IedB,jsd:jed,nz)) ; CS%diffu(:,:,:) = 0.0
1353	1534325	ALLOC_(CS%diffv(isd:ied,JsdB:JedB,nz)) ; CS%diffv(:,:,:) = 0.0
1354		!$omp target enter data map(to: CS%diffu, CS%diffv)
1355	1530343	ALLOC_(CS%CAu(IsdB:IedB,jsd:jed,nz)) ; CS%CAu(:,:,:) = 0.0
1356	1534325	ALLOC_(CS%CAv(isd:ied,JsdB:JedB,nz)) ; CS%CAv(:,:,:) = 0.0
1357		!$omp target enter data map(to: CS%CAu, CS%CAv)
1358	1530343	ALLOC_(CS%CAu_pred(IsdB:IedB,jsd:jed,nz)) ; CS%CAu_pred(:,:,:) = 0.0
1359	1534325	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	1530343	ALLOC_(CS%PFu(IsdB:IedB,jsd:jed,nz)) ; CS%PFu(:,:,:) = 0.0
1362	1534325	ALLOC_(CS%PFv(isd:ied,JsdB:JedB,nz)) ; CS%PFv(:,:,:) = 0.0
1363		!$omp target enter data map(to: CS%PFu, CS%PFv)
1364	69373	ALLOC_(CS%eta(isd:ied,jsd:jed)) ; CS%eta(:,:) = 0.0
1365	1530343	ALLOC_(CS%u_av(IsdB:IedB,jsd:jed,nz)) ; CS%u_av(:,:,:) = 0.0
1366	1534325	ALLOC_(CS%v_av(isd:ied,JsdB:JedB,nz)) ; CS%v_av(:,:,:) = 0.0
1367	1526207	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	69561	allocate(CS%taux_bot(IsdB:IedB,jsd:jed), source=0.0)
1632	69742	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	69561	ALLOC_(CS%uhbt(IsdB:IedB,jsd:jed)) ; CS%uhbt(:,:) = 0.0
1636	69742	ALLOC_(CS%vhbt(isd:ied,JsdB:JedB)) ; CS%vhbt(:,:) = 0.0
1637		!$omp target enter data map(alloc: CS%uhbt, CS%vhbt)
1638	1530343	ALLOC_(CS%visc_rem_u(IsdB:IedB,jsd:jed,nz)) ; CS%visc_rem_u(:,:,:) = 0.0
1639	1534325	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	69373	ALLOC_(CS%eta_PF(isd:ied,jsd:jed)) ; CS%eta_PF(:,:) = 0.0
1642	1526207	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	1530343	ALLOC_(CS%u_accel_bt(IsdB:IedB,jsd:jed,nz)) ; CS%u_accel_bt(:,:,:) = 0.0
1645	1534325	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	1530343	ALLOC_(CS%PFu_Stokes(IsdB:IedB,jsd:jed,nz)) ; CS%PFu_Stokes(:,:,:) = 0.0
1648	1534325	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	64981	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	1429583	do k=1,nz ; do j=js,je ; do i=is,ie
1736	1429560	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	64981	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	1530343	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	1534325	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	1591601	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	1591601	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