← back to index

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.
6module MOM_dynamics_split_RK2
7
8use MOM_variables, only : vertvisc_type, thermo_var_ptrs, porous_barrier_type
9use MOM_variables, only : BT_cont_type, alloc_bt_cont_type, dealloc_bt_cont_type
10use MOM_variables, only : accel_diag_ptrs, ocean_internal_state, cont_diag_ptrs
11use MOM_forcing_type, only : mech_forcing
12
13use MOM_checksum_packages, only : MOM_thermo_chksum, MOM_state_chksum, MOM_accel_chksum
14use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
15use MOM_cpu_clock, only : CLOCK_COMPONENT, CLOCK_SUBCOMPONENT
16use MOM_cpu_clock, only : CLOCK_MODULE_DRIVER, CLOCK_MODULE, CLOCK_ROUTINE
17use MOM_diabatic_driver, only : diabatic_CS, extract_diabatic_member
18use MOM_diag_mediator, only : diag_mediator_init, enable_averages
19use MOM_diag_mediator, only : disable_averaging, post_data, safe_alloc_ptr
20use MOM_diag_mediator, only : post_product_u, post_product_sum_u
21use MOM_diag_mediator, only : post_product_v, post_product_sum_v
22use MOM_diag_mediator, only : register_diag_field, register_static_field
23use MOM_diag_mediator, only : set_diag_mediator_grid, diag_ctrl, diag_update_remap_grids
24use MOM_domains, only : To_South, To_West, To_All, CGRID_NE, SCALAR_PAIR
25use MOM_domains, only : To_North, To_East, Omit_Corners
26use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type
27use MOM_domains, only : start_group_pass, complete_group_pass, pass_var, pass_vector
28use MOM_debugging, only : hchksum, uvchksum, query_debugging_checks
29use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, WARNING, is_root_pe
30use MOM_error_handler, only : MOM_set_verbosity, callTree_showQuery
31use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint
32use MOM_file_parser, only : get_param, log_version, param_file_type
33use MOM_get_input, only : directories
34use MOM_io, only : vardesc, var_desc, EAST_FACE, NORTH_FACE
35use MOM_restart, only : register_restart_field, register_restart_pair
36use MOM_restart, only : query_initialized, set_initialized, save_restart
37use MOM_restart, only : only_read_from_restarts
38use MOM_restart, only : restart_init, is_new_run, MOM_restart_CS
39use MOM_time_manager, only : time_type, operator(+)
40use MOM_time_manager, only : operator(-), operator(>), operator(*), operator(/)
41
42use MOM_ALE, only : ALE_CS, ALE_remap_velocities
43use MOM_barotropic, only : barotropic_init, btstep, btcalc, bt_mass_source
44use MOM_barotropic, only : register_barotropic_restarts, set_dtbt, barotropic_CS
45use MOM_barotropic, only : barotropic_end
46use MOM_boundary_update, only : update_OBC_data, update_OBC_CS
47use MOM_continuity, only : continuity, continuity_CS
48use MOM_continuity, only : continuity_init, continuity_stencil
49use MOM_CoriolisAdv, only : CorAdCalc, CoriolisAdv_CS
50use MOM_CoriolisAdv, only : CoriolisAdv_init, CoriolisAdv_end, CoriolisAdv_stencil
51use MOM_CVMix_KPP, only : KPP_get_BLD, KPP_CS
52use MOM_debugging, only : check_redundant
53use MOM_energetic_PBL, only : energetic_PBL_get_MLD, energetic_PBL_CS
54use MOM_grid, only : ocean_grid_type
55use MOM_harmonic_analysis, only : HA_init, harmonic_analysis_CS
56use MOM_hor_index, only : hor_index_type
57use MOM_hor_visc, only : horizontal_viscosity, hor_visc_CS, hor_visc_vel_stencil
58use MOM_hor_visc, only : hor_visc_init, hor_visc_end
59use MOM_interface_heights, only : thickness_to_dz, find_col_avg_SpV
60use MOM_lateral_mixing_coeffs, only : VarMix_CS
61use MOM_MEKE_types, only : MEKE_type
62use MOM_open_boundary, only : ocean_OBC_type, radiation_open_bdry_conds
63use MOM_open_boundary, only : open_boundary_zero_normal_flow, open_boundary_query
64use MOM_open_boundary, only : open_boundary_test_extern_h, update_OBC_ramp
65use MOM_open_boundary, only : copy_thickness_reservoirs
66use MOM_open_boundary, only : update_segment_thickness_reservoirs
67use MOM_PressureForce, only : PressureForce, PressureForce_CS
68use MOM_PressureForce, only : PressureForce_init
69use MOM_set_visc, only : set_viscous_ML, set_visc_CS
70use MOM_stochastics, only : stochastic_CS
71use MOM_thickness_diffuse, only : thickness_diffuse_CS
72use MOM_self_attr_load, only : SAL_CS
73use MOM_self_attr_load, only : SAL_init, SAL_end
74use MOM_tidal_forcing, only : tidal_forcing_CS
75use MOM_tidal_forcing, only : tidal_forcing_init, tidal_forcing_end
76use MOM_unit_scaling, only : unit_scale_type
77use MOM_vert_friction, only : vertvisc, vertvisc_coef, vertvisc_remnant
78use MOM_vert_friction, only : vertvisc_init, vertvisc_end, vertvisc_CS
79use MOM_vert_friction, only : updateCFLtruncationValue, vertFPmix
80use MOM_verticalGrid, only : verticalGrid_type, get_thickness_units
81use MOM_verticalGrid, only : get_flux_units, get_tr_flux_units
82use MOM_wave_interface, only : wave_parameters_CS, Stokes_PGF
83
84implicit none ; private
85
86#include <MOM_memory.h>
87
88!> MOM_dynamics_split_RK2 module control structure
89type, 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
281end type MOM_dyn_split_RK2_CS
282
283
284public step_MOM_dyn_split_RK2
285public register_restarts_dyn_split_RK2
286public initialize_dyn_split_RK2
287public remap_dyn_split_RK2_aux_vars
288public init_dyn_split_RK2_diabatic
289public end_dyn_split_RK2
290
291!>@{ CPU time clock IDs
292integer :: id_clock_Cor, id_clock_pres, id_clock_vertvisc
293integer :: id_clock_horvisc, id_clock_mom_update
294integer :: id_clock_continuity, id_clock_thick_diff
295integer :: id_clock_btstep, id_clock_btcalc, id_clock_btforce
296integer :: id_clock_pass, id_clock_pass_init
297!>@}
298
299contains
300
301!> RK2 splitting for time stepping MOM adiabatic dynamics
30224subroutine step_MOM_dyn_split_RK2(u_inst, v_inst, h, tv, visc, Time_local, dt, forces, &
30324 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
34948 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: up ! Predicted zonal velocity [L T-1 ~> m s-1].
35048 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: vp ! Predicted meridional velocity [L T-1 ~> m s-1].
35148 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: hp ! Predicted thickness [H ~> m or kg m-2].
35248 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: dz ! Distance between the interfaces around a layer [Z ~> m]
35348 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: ueffA ! Effective Area of U-Faces [H L ~> m2]
35448 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: veffA ! Effective Area of V-Faces [H L ~> m2]
35548 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]
35848 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
36248 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]
36448 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
36748 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]
36948 real, dimension(SZI_(G),SZJ_(G)) :: SpV_avg ! The column averaged specific volume [R-1 ~> m3 kg-1]
37048 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
37348 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]
37548 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?
37948 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]
38148 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.
40224 u_av, & ! The zonal velocity time-averaged over a time step [L T-1 ~> m s-1].
40324 v_av, & ! The meridional velocity time-averaged over a time step [L T-1 ~> m s-1].
40424 h_av ! The layer thickness time-averaged over a time step [H ~> m or kg m-2].
405
40648 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
42148 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
42424 is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
42524 Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
42624 u_av => CS%u_av ; v_av => CS%v_av ; h_av => CS%h_av ; eta => CS%eta
427
42824 Idt_bc = 1.0 / dt
429
43024 sym = G%Domain%symmetric ! switch to include symmetric domain in checksums
431
43224 showCallTree = callTree_showQuery()
4330 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
4383120 do concurrent (k=1:nz, j=G%jsd:G%jed, I=G%IsdB:G%IedB)
43916003248 up(I,j,k) = 0.0
440 enddo
4413096 do concurrent (k=1:nz, J=G%JsdB:G%JedB, i=G%isd:G%ied)
44216112664 vp(i,J,k) = 0.0
443 enddo
4443096 do concurrent (k=1:nz, j=G%jsd:G%jed, i=G%isd:G%ied)
44515879192 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
45024 call updateCFLtruncationValue(Time_local, CS%vertvisc_CSp, US)
451
45224 if (CS%debug) then
4530 call query_debugging_checks(do_redundant=debug_redundant)
4540 call MOM_state_chksum("Start predictor ", u_inst, v_inst, h, uh, vh, G, GV, US, symmetric=sym)
4550 if (debug_redundant) then
4560 call check_redundant("Start predictor u ", u_inst, v_inst, G, unscale=US%L_T_to_m_s)
4570 call check_redundant("Start predictor uh ", uh, vh, G, unscale=GV%H_to_MKS*US%L_to_m**2*US%s_to_T)
458 endif
459 endif
460
46124 dyn_p_surf = associated(p_surf_begin) .and. associated(p_surf_end)
46224 if (dyn_p_surf) then
4630 p_surf => p_surf_end
4640 call safe_alloc_ptr(eta_PF_start,G%isd,G%ied,G%jsd,G%jed)
4650 eta_PF_start(:,:) = 0.0
466 else
46724 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
47424 if (associated(CS%OBC)) then
4750 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
4780 call update_OBC_ramp(Time_local, CS%OBC, US)
479
4800 do k=1,nz ; do j=G%jsd,G%jed ; do I=G%IsdB,G%IedB
4810 u_old_rad_OBC(I,j,k) = u_av(I,j,k)
482 enddo ; enddo ; enddo
4830 do k=1,nz ; do J=G%JsdB,G%JedB ; do i=G%isd,G%ied
4840 v_old_rad_OBC(i,J,k) = v_av(i,J,k)
485 enddo ; enddo ; enddo
486 endif
487
48824 BT_cont_BT_thick = .false.
48924 if (associated(CS%BT_cont)) BT_cont_BT_thick = &
49024 (allocated(CS%BT_cont%h_u) .and. allocated(CS%BT_cont%h_v))
491
49224 if (CS%split_bottom_stress) then
4930 taux_bot => CS%taux_bot ; tauy_bot => CS%tauy_bot
494 endif
495
496 !--- begin set up for group halo pass
497
49824 cont_stencil = continuity_stencil(CS%continuity_CSp)
49924 obc_stencil = 2
50024 cor_stencil = CoriolisAdv_stencil(CS%CoriolisAdv)
50124 if (associated(CS%OBC)) then
5020 if (CS%OBC%oblique_BCs_exist_globally) obc_stencil = 3
503 endif
50424 vel_stencil = max(2, obc_stencil, hor_visc_vel_stencil(CS%hor_visc))
50524 call cpu_clock_begin(id_clock_pass)
50624 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, &
50824 To_All+SCALAR_PAIR, CGRID_NE, halo=max(1,cont_stencil))
50924 call create_group_pass(CS%pass_uvp, up, vp, G%Domain, halo=max(1,cont_stencil))
51024 call create_group_pass(CS%pass_uv, u_inst, v_inst, G%Domain, halo=max(2,cont_stencil))
511
51224 call create_group_pass(CS%pass_hp_uv, hp, G%Domain, halo=cor_stencil)
51324 call create_group_pass(CS%pass_hp_uv, u_av, v_av, G%Domain, halo=max(cor_stencil,vel_stencil))
51424 call create_group_pass(CS%pass_hp_uv, uh(:,:,:), vh(:,:,:), G%Domain, halo=max(cor_stencil,vel_stencil))
51524 call create_group_pass(CS%pass_h, h, g%domain, halo=max(cor_stencil,cont_stencil))
51624 call create_group_pass(CS%pass_av_uvh, u_av, v_av, g%domain, halo=max(cor_stencil,vel_stencil))
51724 call create_group_pass(CS%pass_av_uvh, uh(:,:,:), vh(:,:,:), G%Domain, halo=max(cor_stencil,vel_stencil))
518
51924 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
52424 if (CS%begw == 0.0) call enable_averages(dt, Time_local, CS%diag)
52524 call cpu_clock_begin(id_clock_pres)
526
527 call PressureForce(h, tv, CS%PFu, CS%PFv, G, GV, US, CS%PressureForce_CSp, &
52824 CS%ALE_CSp, CS%ADp, p_surf, CS%pbce, CS%eta_PF)
529
53024 if (dyn_p_surf) then
531 !$omp target update from(CS%eta_PF)
5320 pres_to_eta = 1.0 / (GV%g_Earth * GV%H_to_RZ)
533 !$OMP parallel do default(shared)
5340 do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
5350 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
53924 Use_Stokes_PGF = present(Waves)
54024 if (Use_Stokes_PGF) then
54124 Use_Stokes_PGF = associated(Waves)
54224 if (Use_Stokes_PGF) Use_Stokes_PGF = Waves%Stokes_PGF
54324 if (Use_Stokes_PGF) then
544 !$omp target update from(h, u_inst, v_inst, CS%PFu, CS%PFv)
5450 call thickness_to_dz(h, tv, dz, G, GV, US, halo_size=1)
5460 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.
5510 if (.not.Waves%Passive_Stokes_PGF) then
5520 do k=1,nz
5530 do j=js,je ; do I=Isq,Ieq
5540 CS%PFu(I,j,k) = CS%PFu(I,j,k) + CS%PFu_Stokes(I,j,k)
555 enddo ; enddo
556 enddo
5570 do k=1,nz
5580 do J=Jsq,Jeq ; do i=is,ie
5590 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
56524 call cpu_clock_end(id_clock_pres)
56624 call disable_averaging(CS%diag)
56724 if (showCallTree) call callTree_wayPoint("done with PressureForce (step_MOM_dyn_split_RK2)")
568
56924 if (associated(CS%OBC)) then ; if (CS%OBC%update_OBC) then
5700 call update_OBC_data(CS%OBC, G, GV, US, tv, h, CS%update_OBC_CSp, Time_local)
571 endif ; endif
572
57324 if (associated(CS%OBC) .and. CS%debug_OBC) then
574 !$omp target update from (CS%PFu, CS%PFv)
5750 call open_boundary_zero_normal_flow(CS%OBC, G, GV, CS%PFu, CS%PFv)
576 endif
577
57824 if (G%nonblocking_updates) then
579 !$omp target update from(eta)
5800 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
58424 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.
5870 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, &
5900 G, GV, US, CS%CoriolisAdv, pbv, Waves=Waves)
5910 call cpu_clock_end(id_clock_Cor)
5920 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])
59624 call cpu_clock_begin(id_clock_btforce)
597
59824 do concurrent (k=1:nz, j=js:je, I=Isq:Ieq)
59913245168 u_bc_accel(I,j,k) = (CS%CAu_pred(I,j,k) + CS%PFu(I,j,k)) + CS%diffu(I,j,k)
600 enddo
60124 do concurrent (k=1:nz, J=Jsq:Jeq, i=is:ie)
60213354584 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
60524 if (associated(CS%OBC)) then
606 !$omp target update from(u_bc_accel, v_bc_accel)
6070 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
61024 call cpu_clock_end(id_clock_btforce)
611
61224 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, &
6190 symmetric=sym)
6200 if (debug_redundant) then
6210 call check_redundant("pre-btstep CS%CA ", CS%CAu_pred, CS%CAv_pred, G, unscale=US%L_T2_to_m_s2)
6220 call check_redundant("pre-btstep CS%PF ", CS%PFu, CS%PFv, G, unscale=US%L_T2_to_m_s2)
6230 call check_redundant("pre-btstep CS%diff ", CS%diffu, CS%diffv, G, unscale=US%L_T2_to_m_s2)
6240 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
62824 call cpu_clock_begin(id_clock_vertvisc)
629
63024 do concurrent (k=1:nz, j=js:je, I=Isq:Ieq)
63113245168 up(I,j,k) = G%mask2dCu(I,j) * (u_inst(I,j,k) + dt * u_bc_accel(I,j,k))
632 enddo
63324 do concurrent (k=1:nz, J=Jsq:Jeq, i=is:ie)
63413354584 vp(i,J,k) = G%mask2dCv(i,J) * (v_inst(i,J,k) + dt * v_bc_accel(i,J,k))
635 enddo
636
63724 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.
64024 call set_viscous_ML(u_inst, v_inst, h, tv, forces, visc, dt, G, GV, US, CS%set_visc_CSp)
64124 call disable_averaging(CS%diag)
642
64324 if (CS%debug) then
644 !$omp target update from(up, vp)
6450 call uvchksum("before vertvisc: up", up, vp, G%HI, haloshift=0, symmetric=sym, unscale=US%L_T_to_m_s)
646 endif
647
64824 call thickness_to_dz(h, tv, dz, G, GV, US, halo_size=1, do_offload=.true.)
649
65024 call vertvisc_coef(up, vp, h, dz, forces, visc, tv, dt, G, GV, US, CS%vertvisc_CSp, CS%OBC, VarMix)
65124 call vertvisc_remnant(visc, CS%visc_rem_u, CS%visc_rem_v, dt, G, GV, US, CS%vertvisc_CSp)
652
65324 call cpu_clock_end(id_clock_vertvisc)
65424 if (showCallTree) call callTree_wayPoint("done with vertvisc_coef (step_MOM_dyn_split_RK2)")
655
65624 call cpu_clock_begin(id_clock_pass)
65724 if (G%nonblocking_updates) then
6580 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)
6610 call start_group_pass(CS%pass_visc_rem, G%Domain)
662 else
66324 call do_group_pass(CS%pass_eta, G%Domain, omp_offload=.true.)
66424 call do_group_pass(CS%pass_visc_rem, G%Domain, omp_offload=.true.)
665 endif
66624 call cpu_clock_end(id_clock_pass)
667
66824 call cpu_clock_begin(id_clock_btcalc)
669 ! Calculate the relative layer weights for determining barotropic quantities.
67024 if (.not.BT_cont_BT_thick) then
6710 call btcalc(h, G, GV, CS%barotropic_CSp, OBC=CS%OBC)
672 endif
67324 call bt_mass_source(h, eta, .true., G, GV, CS%barotropic_CSp)
674
675210552 SpV_avg(:,:) = 0.0
67624 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.
6790 if (open_boundary_query(CS%OBC, apply_Flather_OBC=.true.)) &
6800 call find_col_avg_SpV(h, SpV_avg, tv, G, GV, US)
681 endif
68224 call cpu_clock_end(id_clock_btcalc)
683
68424 if (G%nonblocking_updates) then
6850 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
68924 if (associated(CS%OBC)) &
6900 call copy_thickness_reservoirs(CS%OBC, G, GV)
691
692! u_accel_bt = layer accelerations due to barotropic solver
69324 if (associated(CS%BT_cont) .or. CS%BT_use_layer_fluxes) then
69424 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, &
69624 visc_rem_u=CS%visc_rem_u, visc_rem_v=CS%visc_rem_v, BT_cont=CS%BT_cont)
697
69824 call cpu_clock_end(id_clock_continuity)
69924 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, &
70224 OBC=CS%OBC)
703
704 endif
70524 if (showCallTree) call callTree_wayPoint("done with continuity[BT_cont] (step_MOM_dyn_split_RK2)")
706 endif
707
70824 if (CS%BT_use_layer_fluxes) then
70924 uh_ptr => uh_in ; vh_ptr => vh_in ; u_ptr => u_inst ; v_ptr => v_inst
710 endif
711
71224 call cpu_clock_begin(id_clock_btstep)
71324 if (calc_dtbt) then
71413 if (CS%dtbt_use_bt_cont .and. associated(CS%BT_cont)) then
7150 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.
71913 call set_dtbt(G, GV, US, CS%barotropic_CSp, CS%pbce, eta=eta)
720 endif
721 endif
722
72324 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, &
72924 eta_PF_start, taux_bot, tauy_bot, uh_ptr, vh_ptr, u_ptr, v_ptr)
730
73124 if (showCallTree) call callTree_leave("btstep()")
73224 call cpu_clock_end(id_clock_btstep)
733
734! up = u + dt_pred*( u_bc_accel + u_accel_bt )
73524 dt_pred = dt * CS%be
73624 call cpu_clock_begin(id_clock_mom_update)
737
73824 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 * &
74013354584 (v_bc_accel(i,J,k) + CS%v_accel_bt(i,J,k)))
741 enddo
742
74324 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 * &
74513245168 (u_bc_accel(I,j,k) + CS%u_accel_bt(I,j,k)))
746 enddo
747
74824 call cpu_clock_end(id_clock_mom_update)
749
75024 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, &
7550 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)
7580 call uvchksum("Predictor 1 [uv]", up, vp, G%HI, haloshift=0, symmetric=sym, unscale=US%L_T_to_m_s)
7590 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, &
7610 symmetric=sym, unscale=GV%H_to_MKS*US%L_to_m**2*US%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, &
7670 symmetric=sym)
7680 if (debug_redundant) then
7690 call check_redundant("Predictor 1 up", up, vp, G, unscale=US%L_T_to_m_s)
7700 call check_redundant("Predictor 1 uh", uh, vh, G, unscale=GV%H_to_MKS*US%L_to_m**2*US%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
77624 call cpu_clock_begin(id_clock_vertvisc)
77724 if (CS%debug) then
7780 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
78124 if (CS%fpmix) then
782 !$omp target update from(up, vp)
7830 uold(:,:,:) = 0.0
7840 vold(:,:,:) = 0.0
7850 do k = 1, nz
7860 do j = js , je
7870 do I = Isq, Ieq
7880 uold(I,j,k) = up(I,j,k)
789 enddo
790 enddo
7910 do J = Jsq, Jeq
7920 do i = is, ie
7930 vold(i,J,k) = vp(i,J,k)
794 enddo
795 enddo
796 enddo
797 endif
798
79924 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, &
80224 CS%OBC, VarMix)
803
80424 if (CS%fpmix) then
8050 hbl(:,:) = 0.0
8060 if (ASSOCIATED(CS%KPP_CSp)) call KPP_get_BLD(CS%KPP_CSp, hbl, G, US, m_to_BLD_units=GV%m_to_H)
8070 if (ASSOCIATED(CS%energetic_PBL_CSp)) &
8080 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
8110 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, &
8140 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, &
8180 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, &
82124 GV, US, CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot, waves=waves)
822
823 endif
824
82524 if (showCallTree) call callTree_wayPoint("done with vertvisc (step_MOM_dyn_split_RK2)")
82624 if (G%nonblocking_updates) then
8270 call cpu_clock_end(id_clock_vertvisc)
828 !$omp target update from(up, vp)
8290 call start_group_pass(CS%pass_uvp, G%Domain, clock=id_clock_pass)
8300 call cpu_clock_begin(id_clock_vertvisc)
831 endif
832
83324 if (CS%visc_rem_dt_bug) then
8340 call vertvisc_remnant(visc, CS%visc_rem_u, CS%visc_rem_v, dt_pred, G, GV, US, CS%vertvisc_CSp)
835 else
83624 call vertvisc_remnant(visc, CS%visc_rem_u, CS%visc_rem_v, dt, G, GV, US, CS%vertvisc_CSp)
837 endif
838
83924 call cpu_clock_end(id_clock_vertvisc)
840
84124 call do_group_pass(CS%pass_visc_rem, G%Domain, clock=id_clock_pass, omp_offload=.true.)
842
84324 if (G%nonblocking_updates) then
8440 call complete_group_pass(CS%pass_uvp, G%Domain, clock=id_clock_pass)
845 !$omp target update to(up, vp)
846 else
84724 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
85224 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, &
85524 u_cor=u_av, v_cor=v_av, BT_cont=CS%BT_cont)
85624 call cpu_clock_end(id_clock_continuity)
857
85824 if (showCallTree) call callTree_wayPoint("done with continuity (step_MOM_dyn_split_RK2)")
859
86024 call do_group_pass(CS%pass_hp_uv, G%Domain, clock=id_clock_pass, omp_offload=.true.)
861
86224 if (associated(CS%OBC)) then
863 !$omp target update from(u_av, v_av)
864
8650 if (CS%debug) &
8660 call uvchksum("Pre OBC avg [uv]", u_av, v_av, G%HI, haloshift=1, symmetric=sym, unscale=US%L_T_to_m_s)
867
8680 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
8700 if (CS%debug) &
8710 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)
8770 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
8813000 do concurrent (k=1:nz, j=js-cor_stencil:je+cor_stencil, i=is-cor_stencil:ie+cor_stencil)
88214478264 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.
88624 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)
89224 if (CS%BT_adj_corr_mass_src) then
89324 call cpu_clock_begin(id_clock_btcalc)
89424 call bt_mass_source(hp, eta_pred, .false., G, GV, CS%barotropic_CSp)
89524 call cpu_clock_end(id_clock_btcalc)
896 endif
897
89824 if (CS%begw /= 0.0) then
899 ! hp <- (1-begw)*h_in + begw*hp
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.
9020 do concurrent (k=1:nz, j=js-2:je+2, i=is-2:ie+2)
9030 hp(i,j,k) = (1.0-CS%begw)*h(i,j,k) + CS%begw*hp(i,j,k)
904 enddo
905
906 ! PFu = d/dx M(hp,T,S)
907 ! pbce = dM/deta
9080 call cpu_clock_begin(id_clock_pres)
909 call PressureForce(hp, tv, CS%PFu, CS%PFv, G, GV, US, CS%PressureForce_CSp, &
9100 CS%ALE_CSp, CS%ADp, p_surf, CS%pbce, CS%eta_PF)
911 ! Stokes shear force contribution to pressure gradient
9120 Use_Stokes_PGF = present(Waves)
9130 if (Use_Stokes_PGF) then
9140 Use_Stokes_PGF = associated(Waves)
9150 if (Use_Stokes_PGF) Use_Stokes_PGF = Waves%Stokes_PGF
9160 if (Use_Stokes_PGF) then
917 !$omp target update from(CS%PFu, CS%PFv, h)
9180 call thickness_to_dz(h, tv, dz, G, GV, US, halo_size=1)
9190 call Stokes_PGF(G, GV, US, dz, u_inst, v_inst, CS%PFu_Stokes, CS%PFv_Stokes, Waves)
9200 if (.not.Waves%Passive_Stokes_PGF) then
9210 do k=1,nz
9220 do j=js,je ; do I=Isq,Ieq
9230 CS%PFu(I,j,k) = CS%PFu(I,j,k) + CS%PFu_Stokes(I,j,k)
924 enddo ; enddo
925 enddo
9260 do k=1,nz
9270 do J=Jsq,Jeq ; do i=is,ie
9280 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
9340 call cpu_clock_end(id_clock_pres)
9350 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
94124 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, &
94324 OBC=CS%OBC)
94424 if (showCallTree) call callTree_wayPoint("done with btcalc[BT_cont_BT_thick] (step_MOM_dyn_split_RK2)")
945 endif
946
94724 if (CS%debug) then
948 !$omp target update from(up, vp, hp, uh, vh)
9490 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)
9510 call uvchksum("Predictor avg [uv]", u_av, v_av, G%HI, haloshift=1, symmetric=sym, unscale=US%L_T_to_m_s)
9520 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)
9540 if (debug_redundant) then
9550 call check_redundant("Predictor up ", up, vp, G, unscale=US%L_T_to_m_s)
9560 call check_redundant("Predictor uh ", uh, vh, G, unscale=GV%H_to_MKS*US%L_to_m**2*US%s_to_T)
957 endif
958 endif
959
960! diffu = horizontal viscosity terms (u_av)
96124 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, tv, dt, &
964 OBC=CS%OBC, BT=CS%barotropic_CSp, TD=thickness_diffuse_CSp, &
96524 ADp=CS%ADp, hu_cont=CS%BT_cont%h_u, hv_cont=CS%BT_cont%h_v, STOCH=STOCH)
966
96724 call cpu_clock_end(id_clock_horvisc)
96824 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
97124 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, &
97324 G, GV, US, CS%CoriolisAdv, pbv, Waves=Waves)
97424 call cpu_clock_end(id_clock_Cor)
975
97624 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])
98124 call cpu_clock_begin(id_clock_btforce)
982
98324 do concurrent (k=1:nz, j=js:je, I=Isq:Ieq)
98413245168 u_bc_accel(I,j,k) = (CS%Cau(I,j,k) + CS%PFu(I,j,k)) + CS%diffu(I,j,k)
985 enddo
986
98724 do concurrent (k=1:nz, J=Jsq:Jeq,i=is:ie)
98813354584 v_bc_accel(i,J,k) = (CS%Cav(i,J,k) + CS%PFv(i,J,k)) + CS%diffv(i,J,k)
989 enddo
990
99124 if (associated(CS%OBC)) then
992 !$omp target update from(u_bc_accel, v_bc_accel)
9930 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
99624 call cpu_clock_end(id_clock_btforce)
997
99824 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, &
10050 symmetric=sym)
10060 if (debug_redundant) then
10070 call check_redundant("corr pre-btstep CS%CA ", CS%CAu, CS%CAv, G, unscale=US%L_T2_to_m_s2)
10080 call check_redundant("corr pre-btstep CS%PF ", CS%PFu, CS%PFv, G, unscale=US%L_T2_to_m_s2)
10090 call check_redundant("corr pre-btstep CS%diff ", CS%diffu, CS%diffv, G, unscale=US%L_T2_to_m_s2)
10100 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
101624 call cpu_clock_begin(id_clock_btstep)
101724 if (CS%BT_use_layer_fluxes) then
101824 uh_ptr => uh ; vh_ptr => vh ; u_ptr => u_av ; v_ptr => v_av
1019 endif
1020
102124 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, &
102624 eta_PF_start, taux_bot, tauy_bot, uh_ptr, vh_ptr, u_ptr, v_ptr, etaav=eta_av)
1027
102824 if (CS%id_deta_dt>0) then
1029 !$omp target update from(eta, eta_pred)
10300 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
103324 do concurrent (j=js:je, i=is:ie)
1034175704 eta(i,j) = eta_pred(i,j)
1035 enddo
1036
103724 call cpu_clock_end(id_clock_btstep)
103824 if (showCallTree) call callTree_leave("btstep()")
1039
104024 if (CS%debug .and. debug_redundant) then
1041 !$omp target update from(CS%u_accel_bt, CS%v_accel_bt)
10420 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 )
104624 call cpu_clock_begin(id_clock_mom_update)
104724 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 * &
104913245168 (u_bc_accel(I,j,k) + CS%u_accel_bt(I,j,k)))
1050 enddo
105124 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 * &
105313354584 (v_bc_accel(i,J,k) + CS%v_accel_bt(i,J,k)))
1054 enddo
105524 call cpu_clock_end(id_clock_mom_update)
1056
105724 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, &
10640 symmetric=sym)
1065
1066 !$omp target update from(u_inst, v_inst, h, uh, vh)
10670 call uvchksum("Corrector 1 [uv]", u_inst, v_inst, G%HI, haloshift=0, symmetric=sym, unscale=US%L_T_to_m_s)
10680 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, &
10700 symmetric=sym, unscale=GV%H_to_MKS*US%L_to_m**2*US%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
107724 call cpu_clock_begin(id_clock_vertvisc)
1078
107924 if (CS%fpmix) then
1080 !$omp target update from(u_inst, v_inst)
10810 uold(:,:,:) = 0.0
10820 vold(:,:,:) = 0.0
10830 do k = 1, nz
10840 do j = js , je
10850 do I = Isq, Ieq
10860 uold(I,j,k) = u_inst(I,j,k)
1087 enddo
1088 enddo
10890 do J = Jsq, Jeq
10900 do i = is, ie
10910 vold(i,J,k) = v_inst(i,J,k)
1092 enddo
1093 enddo
1094 enddo
1095 endif
1096
109724 call thickness_to_dz(h, tv, dz, G, GV, US, halo_size=1, do_offload=.true.)
1098
109924 call vertvisc_coef(u_inst, v_inst, h, dz, forces, visc, tv, dt, G, GV, US, CS%vertvisc_CSp, CS%OBC, VarMix)
1100
110124 if (CS%fpmix) then
11020 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, &
11050 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, &
11090 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, &
111224 CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot, waves=waves)
1113 endif
1114
111524 if (G%nonblocking_updates) then
1116 !$omp target update from(u_inst, v_inst)
11170 call cpu_clock_end(id_clock_vertvisc)
11180 call start_group_pass(CS%pass_uv, G%Domain, clock=id_clock_pass)
11190 call cpu_clock_begin(id_clock_vertvisc)
1120 endif
112124 call vertvisc_remnant(visc, CS%visc_rem_u, CS%visc_rem_v, dt, G, GV, US, CS%vertvisc_CSp)
1122
112324 call cpu_clock_end(id_clock_vertvisc)
112424 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.
11273000 do concurrent (k=1:nz, j=js-cor_stencil:je+cor_stencil, i=is-cor_stencil:ie+cor_stencil)
112814478264 h_av(i,j,k) = h(i,j,k)
1129 enddo
1130
113124 call do_group_pass(CS%pass_visc_rem, G%Domain, clock=id_clock_pass, omp_offload=.true.)
1132
113324 if (G%nonblocking_updates) then
11340 call complete_group_pass(CS%pass_uv, G%Domain, clock=id_clock_pass)
1135 !$omp target update to(u_inst, v_inst)
1136 else
113724 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.
114324 call cpu_clock_begin(id_clock_continuity)
11443096 do concurrent (k=1:nz, j=G%jsd:G%jed, i=G%isd:G%ied)
114515879192 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, &
115024 u_cor=u_av, v_cor=v_av)
115124 call cpu_clock_end(id_clock_continuity)
1152
115324 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.
115824 call diag_update_remap_grids(CS%diag)
115924 if (showCallTree) call callTree_wayPoint("done with continuity (step_MOM_dyn_split_RK2)")
1160
116124 if (G%nonblocking_updates) then
1162 !$omp target update from(u_av, v_av, uh, vh)
11630 call start_group_pass(CS%pass_av_uvh, G%Domain, clock=id_clock_pass)
1164 else
116524 call do_group_pass(CS%pass_av_uvh, G%domain, clock=id_clock_pass, omp_offload=.true.)
1166 endif
1167
116824 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)
11720 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.
11763000 do concurrent (k=1:nz, j=js-cor_stencil:je+cor_stencil, i=is-cor_stencil:ie+cor_stencil)
117714478264 h_av(i,j,k) = 0.5*(h_av(i,j,k) + h(i,j,k))
1178 enddo
1179
118024 if (G%nonblocking_updates) then
11810 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
11853024 do concurrent (k=1:nz, j=js-cor_stencil:je+cor_stencil, I=Isq-cor_stencil:Ieq+cor_stencil)
118614595024 uhtr(I,j,k) = uhtr(I,j,k) + uh(I,j,k) * dt
1187 enddo
1188
11893000 do concurrent (k=1:nz, J=Jsq-cor_stencil:Jeq+cor_stencil, i=is-cor_stencil:ie+cor_stencil)
119014704440 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
119624 if (associated(CS%OBC)) then
11970 call update_segment_thickness_reservoirs(G, GV, uhtr, vhtr, h, CS%OBC)
1198 endif
1199
120024 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.
120324 call cpu_clock_begin(id_clock_Cor)
120424 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, &
120724 G, GV, US, CS%CoriolisAdv, pbv, Waves=Waves)
120824 CS%CAu_pred_stored = .true.
120924 call enable_averages(dt, Time_local, CS%diag) ! Reenable the averaging
121024 call cpu_clock_end(id_clock_Cor)
121124 if (showCallTree) call callTree_wayPoint("done with CorAdCalc (step_MOM_dyn_split_RK2)")
1212 else
12130 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
121924 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.
122324 if (CS%id_PFu > 0) call post_data(CS%id_PFu, CS%PFu, CS%diag)
122424 if (CS%id_PFv > 0) call post_data(CS%id_PFv, CS%PFv, CS%diag)
122524 if (CS%id_CAu > 0) call post_data(CS%id_CAu, CS%CAu, CS%diag)
122624 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.
122924 if (CS%id_uh > 0) call post_data(CS%id_uh, uh, CS%diag)
123024 if (CS%id_vh > 0) call post_data(CS%id_vh, vh, CS%diag)
123124 if (CS%id_uav > 0) call post_data(CS%id_uav, u_av, CS%diag)
123224 if (CS%id_vav > 0) call post_data(CS%id_vav, v_av, CS%diag)
123324 if (CS%id_u_BT_accel > 0) call post_data(CS%id_u_BT_accel, CS%u_accel_bt, CS%diag)
123424 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
123724 if (CS%id_ueffA > 0) then
12380 ueffA(:,:,:) = 0
12390 do k=1,nz ; do j=js,je ; do I=Isq,Ieq
12400 if (abs(up(I,j,k)) > 0.) ueffA(I,j,k) = uh(I,j,k) / up(I,j,k)
1241 enddo ; enddo ; enddo
12420 call post_data(CS%id_ueffA, ueffA, CS%diag)
1243 endif
1244
124524 if (CS%id_veffA > 0) then
12460 veffA(:,:,:) = 0
12470 do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
12480 if (abs(vp(i,J,k)) > 0.) veffA(i,J,k) = vh(i,J,k) / vp(i,J,k)
1249 enddo ; enddo ; enddo
12500 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
126824 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)
126924 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
127224 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)
127324 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
127624 if (CS%id_h_PFu > 0) call post_product_u(CS%id_h_PFu, CS%PFu, CS%ADp%diag_hu, G, nz, CS%diag)
127724 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
128024 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)
128124 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)
128224 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)
128324 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)
128424 if (CS%id_h_CAu > 0) call post_product_u(CS%id_h_CAu, CS%CAu, CS%ADp%diag_hu, G, nz, CS%diag)
128524 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
128824 if (CS%id_intz_u_BT_accel_2d > 0) &
12890 call post_product_sum_u(CS%id_intz_u_BT_accel_2d, CS%u_accel_bt, CS%ADp%diag_hu, G, nz, CS%diag)
129024 if (CS%id_intz_v_BT_accel_2d > 0) &
12910 call post_product_sum_v(CS%id_intz_v_BT_accel_2d, CS%v_accel_bt, CS%ADp%diag_hv, G, nz, CS%diag)
129224 if (CS%id_hf_u_BT_accel_2d > 0) &
12930 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)
129424 if (CS%id_hf_v_BT_accel_2d > 0) &
12950 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)
129624 if (CS%id_h_u_BT_accel > 0) &
12970 call post_product_u(CS%id_h_u_BT_accel, CS%u_accel_bt, CS%ADp%diag_hu, G, nz, CS%diag)
129824 if (CS%id_h_v_BT_accel > 0) &
12990 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],
130224 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)
130324 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)
130424 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)
130524 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)
130624 if (CS%id_u_BT_accel_visc_rem > 0) &
13070 call post_product_u(CS%id_u_BT_accel_visc_rem, CS%u_accel_bt, CS%ADp%visc_rem_u, G, nz, CS%diag)
130824 if (CS%id_v_BT_accel_visc_rem > 0) &
13090 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
131224 if (CS%id_deta_dt > 0) call post_data(CS%id_deta_dt, deta_dt, CS%diag)
1313
131424 if (CS%debug) then
1315 !$omp target update from(u_inst, v_inst, h, uh, vh, u_av, v_av, h_av)
13160 call MOM_state_chksum("Corrector ", u_inst, v_inst, h, uh, vh, G, GV, US, symmetric=sym)
13170 call uvchksum("Corrector avg [uv]", u_av, v_av, G%HI, haloshift=1, symmetric=sym, unscale=US%L_T_to_m_s)
13180 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
132224 if (showCallTree) call callTree_leave("step_MOM_dyn_split_RK2()")
1323
132424end 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.
13291subroutine 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.
134241 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
13471 isd = HI%isd ; ied = HI%ied ; jsd = HI%jsd ; jed = HI%jed ; nz = GV%ke
13481 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:)
1352663076 ALLOC_(CS%diffu(IsdB:IedB,jsd:jed,nz)) ; CS%diffu(:,:,:) = 0.0
1353667651 ALLOC_(CS%diffv(isd:ied,JsdB:JedB,nz)) ; CS%diffv(:,:,:) = 0.0
1354 !$omp target enter data map(to: CS%diffu, CS%diffv)
1355663076 ALLOC_(CS%CAu(IsdB:IedB,jsd:jed,nz)) ; CS%CAu(:,:,:) = 0.0
1356667651 ALLOC_(CS%CAv(isd:ied,JsdB:JedB,nz)) ; CS%CAv(:,:,:) = 0.0
1357 !$omp target enter data map(to: CS%CAu, CS%CAv)
1358663076 ALLOC_(CS%CAu_pred(IsdB:IedB,jsd:jed,nz)) ; CS%CAu_pred(:,:,:) = 0.0
1359667651 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)
1361663076 ALLOC_(CS%PFu(IsdB:IedB,jsd:jed,nz)) ; CS%PFu(:,:,:) = 0.0
1362667651 ALLOC_(CS%PFv(isd:ied,JsdB:JedB,nz)) ; CS%PFv(:,:,:) = 0.0
1363 !$omp target enter data map(to: CS%PFu, CS%PFv)
13648773 ALLOC_(CS%eta(isd:ied,jsd:jed)) ; CS%eta(:,:) = 0.0
1365663076 ALLOC_(CS%u_av(IsdB:IedB,jsd:jed,nz)) ; CS%u_av(:,:,:) = 0.0
1366667651 ALLOC_(CS%v_av(isd:ied,JsdB:JedB,nz)) ; CS%v_av(:,:,:) = 0.0
1367657976 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
13701 thickness_units = get_thickness_units(GV)
13711 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.", &
13761 default=.true., do_not_log=.true.)
1377
13781 if (GV%Boussinesq) then
1379 call register_restart_field(CS%eta, "sfc", .false., restart_CS, &
13801 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, &
13830 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.
13886 vd(1) = var_desc("u2", "m s-1", "Auxiliary Zonal velocity", 'u', 'L')
13896 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, &
13911 conversion=US%L_T_to_m_s)
1392
13931 if (CS%store_CAu) then
13946 vd(1) = var_desc("CAu", "m s-2", "Zonal Coriolis and advactive acceleration", 'u', 'L')
13956 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, &
13971 conversion=US%L_T2_to_m_s2)
1398 else
1399 call register_restart_field(CS%h_av, "h2", .false., restart_CS, &
14000 longname="Auxiliary Layer Thickness", units=thickness_units, conversion=GV%H_to_mks)
1401
14020 vd(1) = var_desc("uh", flux_units, "Zonal thickness flux", 'u', 'L')
14030 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, &
14050 conversion=GV%H_to_MKS*US%L_to_m**2*US%s_to_T)
1406 endif
1407
14086 vd(1) = var_desc("diffu", "m s-2", "Zonal horizontal viscous acceleration", 'u', 'L')
14096 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, &
14111 conversion=US%L_T2_to_m_s2)
1412
14131 call register_barotropic_restarts(HI, GV, US, param_file, CS%barotropic_CSp, restart_CS)
1414
141522end 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.
14190subroutine 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
14370 if (.not.CS%remap_aux) return
1438
14390 if (CS%store_CAu) then
14400 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)
14410 call pass_vector(CS%u_av, CS%v_av, G%Domain, complete=.false.)
14420 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)
14430 call pass_vector(CS%CAu_pred, CS%CAv_pred, G%Domain, complete=.true.)
1444 endif
1445
14460 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
1448end 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.
14520subroutine 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
14560 call extract_diabatic_member(diabatic_CSp, KPP_CSp=CS%KPP_CSp)
14570 call extract_diabatic_member(diabatic_CSp, energetic_PBL_CSp=CS%energetic_PBL_CSp)
1458
14590end 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.
14631subroutine 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
15181 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
15341 is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
15351 isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
15361 IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB
1537
15381 if (.not.associated(CS)) call MOM_error(FATAL, &
15390 "initialize_dyn_split_RK2 called with an unassociated control structure.")
15401 if (CS%module_is_initialized) then
1541 call MOM_error(WARNING, "initialize_dyn_split_RK2 called with a control "// &
15420 "structure that has already been initialized.")
15430 return
1544 endif
15451 CS%module_is_initialized = .true.
1546
15471 CS%diag => diag
1548
15491 call log_version(param_file, mdl, version, "")
1550 call get_param(param_file, mdl, "TIDES", CS%use_tides, &
15511 "If true, apply tidal momentum forcing.", default=.false.)
1552 call get_param(param_file, mdl, "CALCULATE_SAL", CS%calculate_SAL, &
15531 "If true, calculate self-attraction and loading.", default=CS%use_tides)
1554 call get_param(param_file, mdl, "USE_HA", CS%use_HA, &
15551 "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.", &
15581 default=0, do_not_log=.not.CS%use_HA)
15591 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 "//&
15671 "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.", &
15751 units="nondim", default=0.0)
1576 call get_param(param_file, mdl, "SET_DTBT_USE_BT_CONT", CS%dtbt_use_bt_cont, &
15771 "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 "//&
15801 "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 "//&
15841 "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, "//&
15891 "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.", &
15931 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.", &
15961 default=.false.)
15971 if (CS%fpmix) then
1598 call get_param(param_file, "MOM", "CEMP_NL", CS%Cemp_NL, &
1599 "Empirical coefficient of non-local momentum mixing.", &
16000 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. "//&
16061 "The default should be changed to true.", default=.false., do_not_log=.true.)
16071 if (CS%remap_aux .and. .not.CS%store_CAu) call MOM_error(FATAL, &
16080 "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.", &
16111 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.", &
16151 default=.false., old_name="DEBUG_OBC", debuggingParam=.true., do_not_log=.true.)
1616 call get_param(param_file, mdl, "DEBUG_TRUNCATIONS", debug_truncations, &
16171 default=.false.)
1618 call get_param(param_file, mdl, "ENABLE_BUGS_BY_DEFAULT", enable_bugs, &
16191 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 "//&
16241 "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 "//&
16291 "is set by VISC_REM_BUG", default=visc_rem_bug)
1630
16318841 allocate(CS%taux_bot(IsdB:IedB,jsd:jed), source=0.0)
16328902 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
16358841 ALLOC_(CS%uhbt(IsdB:IedB,jsd:jed)) ; CS%uhbt(:,:) = 0.0
16368902 ALLOC_(CS%vhbt(isd:ied,JsdB:JedB)) ; CS%vhbt(:,:) = 0.0
1637 !$omp target enter data map(alloc: CS%uhbt, CS%vhbt)
1638663076 ALLOC_(CS%visc_rem_u(IsdB:IedB,jsd:jed,nz)) ; CS%visc_rem_u(:,:,:) = 0.0
1639667651 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)
16418773 ALLOC_(CS%eta_PF(isd:ied,jsd:jed)) ; CS%eta_PF(:,:) = 0.0
1642657976 ALLOC_(CS%pbce(isd:ied,jsd:jed,nz)) ; CS%pbce(:,:,:) = 0.0
1643 !$omp target enter data map(to: CS%pbce, CS%eta_PF)
1644663076 ALLOC_(CS%u_accel_bt(IsdB:IedB,jsd:jed,nz)) ; CS%u_accel_bt(:,:,:) = 0.0
1645667651 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)
1647663076 ALLOC_(CS%PFu_Stokes(IsdB:IedB,jsd:jed,nz)) ; CS%PFu_Stokes(:,:,:) = 0.0
1648667651 ALLOC_(CS%PFv_Stokes(isd:ied,JsdB:JedB,nz)) ; CS%PFv_Stokes(:,:,:) = 0.0
1649
16501 MIS%diffu => CS%diffu
16511 MIS%diffv => CS%diffv
16521 MIS%PFu => CS%PFu
16531 MIS%PFv => CS%PFv
16541 MIS%CAu => CS%CAu
16551 MIS%CAv => CS%CAv
16561 MIS%pbce => CS%pbce
16571 MIS%u_accel_bt => CS%u_accel_bt
16581 MIS%v_accel_bt => CS%v_accel_bt
16591 MIS%u_av => CS%u_av
16601 MIS%v_av => CS%v_av
1661
16621 CS%ADp => Accel_diag
16631 CS%CDp => Cont_diag
16641 Accel_diag%diffu => CS%diffu
16651 Accel_diag%diffv => CS%diffv
16661 Accel_diag%PFu => CS%PFu
16671 Accel_diag%PFv => CS%PFv
16681 Accel_diag%CAu => CS%CAu
16691 Accel_diag%CAv => CS%CAv
16701 Accel_diag%u_accel_bt => CS%u_accel_bt
16711 Accel_diag%v_accel_bt => CS%v_accel_bt
1672
16731 allocate(CS%AD_pred)
16741 CS%AD_pred%diffu => CS%diffu
16751 CS%AD_pred%diffv => CS%diffv
16761 CS%AD_pred%PFu => CS%PFu
16771 CS%AD_pred%PFv => CS%PFv
16781 CS%AD_pred%CAu => CS%CAu_pred
16791 CS%AD_pred%CAv => CS%CAv_pred
16801 CS%AD_pred%u_accel_bt => CS%u_accel_bt
16811 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
16871 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)
16901 call continuity_init(Time, G, GV, US, param_file, diag, CS%continuity_CSp, CS%OBC)
16911 cont_stencil = continuity_stencil(CS%continuity_CSp)
16921 call CoriolisAdv_init(Time, G, GV, US, param_file, diag, CS%ADp, CS%CoriolisAdv)
16931 cor_stencil = CoriolisAdv_stencil(CS%CoriolisAdv)
16941 dyn_h_stencil = max(cont_stencil, CoriolisAdv_stencil(CS%CoriolisAdv))
16951 if (CS%calculate_SAL) call SAL_init(h, tv, G, GV, US, param_file, CS%SAL_CSp, restart_CS)
16961 if (CS%use_tides) call tidal_forcing_init(Time, G, US, param_file, CS%tides_CSp)
16971 if (CS%use_HA) then
16980 call HA_init(Time, US, param_file, nc, CS%HA_CSp)
16990 HA_CSp => CS%HA_CSp
1700 else
17011 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, &
17061 CS%SAL_CSp, CS%tides_CSp)
1707
1708 !$omp target enter data map(alloc: CS%hor_visc)
17091 call hor_visc_init(Time, G, GV, US, param_file, diag, CS%hor_visc, ADp=CS%ADp)
1710
17111 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, &
17141 ntrunc, CS%vertvisc_CSp, CS%fpmix)
17151 CS%set_visc_CSp => set_visc
17161 call updateCFLtruncationValue(Time, CS%vertvisc_CSp, US, activate=is_new_run(restart_CS) )
1717
17181 if (associated(ALE_CSp)) CS%ALE_CSp => ALE_CSp
17191 if (associated(OBC)) then
17200 CS%OBC => OBC
17210 if (OBC%ramp) call update_OBC_ramp(Time, CS%OBC, US, activate=is_new_run(restart_CS) )
1722 endif
17231 if (associated(update_OBC_CSp)) CS%update_OBC_CSp => update_OBC_CSp
1724
17251 eta_rest_name = "sfc" ; if (.not.GV%Boussinesq) eta_rest_name = "p_bot"
17261 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.
17321 if (GV%Boussinesq) then
17337261 do j=js,je ; do i=is,ie ; CS%eta(i,j) = -GV%Z_to_H * G%bathyT(i,j) ; enddo ; enddo
1734 endif
1735544576 do k=1,nz ; do j=js,je ; do i=is,ie
1736544500 CS%eta(i,j) = CS%eta(i,j) + h(i,j,k)
1737 enddo ; enddo ; enddo
17381 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.
17427261 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, &
17471 CS%OBC, CS%SAL_CSp, HA_CSp)
1748
17491 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 tv, dt, OBC=CS%OBC, BT=CS%barotropic_CSp, TD=thickness_diffuse_CSp, &
17541 hu_cont=CS%BT_cont%h_u, hv_cont=CS%BT_cont%h_v)
1755 !$omp target update from(CS%diffu, CS%diffv)
17561 call set_initialized(CS%diffu, "diffu", restart_CS)
17571 call set_initialized(CS%diffv, "diffv", restart_CS)
1758 endif
1759
17601 if (.not. query_initialized(CS%u_av, "u2", restart_CS) .or. &
1761 .not. query_initialized(CS%v_av, "v2", restart_CS)) then
1762663076 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
1763667651 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
17641 call set_initialized(CS%u_av, "u2", restart_CS)
17651 call set_initialized(CS%v_av, "v2", restart_CS)
1766 endif
1767
1768 !$omp target enter data map(alloc: h_tmp )
17691 if (CS%store_CAu) then
17701 if (query_initialized(CS%CAu_pred, "CAu", restart_CS) .and. &
1771 query_initialized(CS%CAv_pred, "CAv", restart_CS)) then
17720 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, &
17761 success=read_uv, scale=US%m_to_L**2*US%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, &
17791 success=read_h2, scale=1.0/GV%H_to_mks)
17801 if (read_uv .and. read_h2) then
17810 call pass_var(CS%h_av, G%Domain, clock=id_clock_pass_init)
1782 else
17831 do concurrent (k=1:nz, j=jsd:jed, i=isd:ied)
1784661633 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)
17871 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)
17911 call pass_var(h_tmp, G%Domain, clock=id_clock_pass_init)
1792 !$omp target update to(h_tmp)
1793
17941 do concurrent (k=1:nz, j=jsd:jed, i=isd:ied)
1795661633 CS%h_av(i,j,k) = 0.5*(h(i,j,k) + h_tmp(i,j,k))
1796 enddo
1797 endif
17981 call pass_vector(CS%u_av, CS%v_av, G%Domain, halo=cor_stencil, clock=id_clock_pass_init, complete=.false.)
17991 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, &
18021 G, GV, US, CS%CoriolisAdv, pbv) !, Waves=Waves)
1803 !$omp target update from(CS%CAu_pred, CS%CAv_pred)
18041 CS%CAu_pred_stored = .true.
1805 endif
1806 else
18070 CS%CAu_pred_stored = .false.
1808 ! This call is just here to initialize uh and vh.
18090 if (.not. query_initialized(uh, "uh", restart_CS) .or. &
1810 .not. query_initialized(vh, "vh", restart_CS)) then
18110 do k=1,nz ; do j=jsd,jed ; do i=isd,ied ; h_tmp(i,j,k) = h(i,j,k) ; enddo ; enddo ; enddo
18120 call continuity(u, v, h, h_tmp, uh, vh, dt, G, GV, US, CS%continuity_CSp, CS%OBC, pbv)
18130 call pass_var(h_tmp, G%Domain, clock=id_clock_pass_init)
18140 do k=1,nz ; do j=jsd,jed ; do i=isd,ied
18150 CS%h_av(i,j,k) = 0.5*(h(i,j,k) + h_tmp(i,j,k))
1816 enddo ; enddo ; enddo
18170 call set_initialized(uh, "uh", restart_CS)
18180 call set_initialized(vh, "vh", restart_CS)
18190 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, &
18240 success=read_uv, scale=US%m_s_to_L_T*US%T_to_s)
18250 CS%CAu_pred_stored = read_uv
1826 else
18270 if (.not. query_initialized(CS%h_av, "h2", restart_CS)) then
18280 CS%h_av(:,:,:) = h(:,:,:)
18290 call set_initialized(CS%h_av, "h2", restart_CS)
1830 endif
1831 endif
1832 endif
1833 !$omp target exit data map(delete: h_tmp )
18341 call cpu_clock_begin(id_clock_pass_init)
18351 call create_group_pass(pass_av_h_uvh, CS%u_av, CS%v_av, G%Domain, halo=2)
18361 if (CS%CAu_pred_stored) then
18371 call create_group_pass(pass_av_h_uvh, CS%CAu_pred, CS%CAv_pred, G%Domain, halo=2)
1838 else
18390 call create_group_pass(pass_av_h_uvh, CS%h_av, G%Domain, halo=2)
18400 call create_group_pass(pass_av_h_uvh, uh, vh, G%Domain, halo=2)
1841 endif
18421 call do_group_pass(pass_av_h_uvh, G%Domain)
18431 call cpu_clock_end(id_clock_pass_init)
1844
18451 flux_units = get_flux_units(GV)
18461 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_MKS*US%L_to_m**2*US%s_to_T, &
18491 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_MKS*US%L_to_m**2*US%s_to_T, &
18521 x_cell_method='sum', v_extensive=.true.)
1853
1854 CS%id_CAu = register_diag_field('ocean_model', 'CAu', diag%axesCuL, Time, &
18551 '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, &
18571 '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, &
18591 '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, &
18611 '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, &
18641 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, &
18671 x_cell_method='sum', v_extensive=.true.)
18681 if (GV%Boussinesq) then
1869 CS%id_deta_dt = register_diag_field('ocean_model', 'deta_dt', diag%axesT1, Time, &
18701 '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', &
18740 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', &
18991 'm s-2', conversion=US%L_T2_to_m_s2)
19001 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', &
19041 'm s-2', conversion=US%L_T2_to_m_s2)
19051 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', &
19091 'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
19101 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', &
19141 'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
19151 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', &
19191 'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
19201 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', &
19241 'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
19251 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', &
19291 'm s-2', conversion=US%L_T2_to_m_s2)
19301 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', &
19341 'm s-2', conversion=US%L_T2_to_m_s2)
19351 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', &
19391 'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
19401 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', &
19441 'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
19451 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', &
19491 'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
19501 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', &
19541 'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
19551 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, &
19581 '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, &
19601 '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, &
19631 '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, &
19651 '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', &
19791 'm s-2', conversion=US%L_T2_to_m_s2)
19801 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', &
19841 'm s-2', conversion=US%L_T2_to_m_s2)
19851 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', &
19891 'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
19901 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', &
19941 'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
19951 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', &
19991 'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
20001 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', &
20041 'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
20051 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', &
20091 'm s-2', conversion=US%L_T2_to_m_s2)
20101 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', &
20131 'm s-2', conversion=US%L_T2_to_m_s2)
20141 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', &
20181 'm s-2', conversion=US%L_T2_to_m_s2)
20191 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', &
20221 'm s-2', conversion=US%L_T2_to_m_s2)
20231 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', &
20271 'm s-2', conversion=US%L_T2_to_m_s2)
20281 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', &
20311 'm s-2', conversion=US%L_T2_to_m_s2)
20321 if (CS%id_v_BT_accel_visc_rem > 0) call safe_alloc_ptr(CS%ADp%visc_rem_v,isd,ied,JsdB,JedB,nz)
2033
20341 id_clock_Cor = cpu_clock_id('(Ocean Coriolis & mom advection)', grain=CLOCK_MODULE)
20351 id_clock_continuity = cpu_clock_id('(Ocean continuity equation)', grain=CLOCK_MODULE)
20361 id_clock_pres = cpu_clock_id('(Ocean pressure force)', grain=CLOCK_MODULE)
20371 id_clock_vertvisc = cpu_clock_id('(Ocean vertical viscosity)', grain=CLOCK_MODULE)
20381 id_clock_horvisc = cpu_clock_id('(Ocean horizontal viscosity)', grain=CLOCK_MODULE)
20391 id_clock_mom_update = cpu_clock_id('(Ocean momentum increments)', grain=CLOCK_MODULE)
20401 id_clock_pass = cpu_clock_id('(Ocean message passing)', grain=CLOCK_MODULE)
20411 id_clock_btcalc = cpu_clock_id('(Ocean barotropic mode calc)', grain=CLOCK_MODULE)
20421 id_clock_btstep = cpu_clock_id('(Ocean barotropic mode stepping)', grain=CLOCK_MODULE)
20431 id_clock_btforce = cpu_clock_id('(Ocean barotropic forcing calc)', grain=CLOCK_MODULE)
2044
20452end subroutine initialize_dyn_split_RK2
2046
2047
2048!> Close the dyn_split_RK2 module
20491subroutine 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)
20531 call barotropic_end(CS%barotropic_CSp)
2054
20551 call vertvisc_end(CS%vertvisc_CSp)
20561 deallocate(CS%vertvisc_CSp)
2057
20581 call hor_visc_end(CS%hor_visc)
2059 !$omp target exit data map(delete: CS%hor_visc)
2060
20611 if (CS%calculate_SAL) call SAL_end(CS%SAL_CSp)
20621 if (CS%use_tides) call tidal_forcing_end(CS%tides_CSp)
20631 call CoriolisAdv_end(CS%CoriolisAdv)
2064
20651 DEALLOC_(CS%diffu) ; DEALLOC_(CS%diffv)
2066 !$omp target exit data map(delete: CS%diffu, CS%diffv)
20671 DEALLOC_(CS%CAu) ; DEALLOC_(CS%CAv)
2068 !$omp target exit data map(delete: CS%CAu, CS%CAv)
20691 DEALLOC_(CS%CAu_pred) ; DEALLOC_(CS%CAv_pred)
2070 !$omp target exit data map(delete: CS%CAu_pred, CS%CAv_pred)
20711 DEALLOC_(CS%PFu) ; DEALLOC_(CS%PFv)
2072 !$omp target exit data map(delete: CS%PFu, CS%PFv)
2073
20741 if (associated(CS%taux_bot)) deallocate(CS%taux_bot)
20751 if (associated(CS%tauy_bot)) deallocate(CS%tauy_bot)
20761 DEALLOC_(CS%uhbt) ; DEALLOC_(CS%vhbt)
20771 DEALLOC_(CS%u_accel_bt) ; DEALLOC_(CS%v_accel_bt)
20781 DEALLOC_(CS%visc_rem_u) ; DEALLOC_(CS%visc_rem_v)
2079
20801 DEALLOC_(CS%eta) ; DEALLOC_(CS%eta_PF) ; DEALLOC_(CS%pbce)
2081 !$omp target exit data map(delete: CS%eta, CS%eta_PF, CS%pbce)
20821 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
20851 call dealloc_BT_cont_type(CS%BT_cont)
20861 deallocate(CS%AD_pred)
2087
20881 deallocate(CS)
20891end 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
21180end module MOM_dynamics_split_RK2