MOM_dynamics_split_RK2b.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> Time step the adiabatic dynamic core of MOM using RK2 method with greater use of the
6!! time-filtered velocities and less inheritance of tedencies from the previous step in the
7!! predictor step than in the original MOM_dyanmics_split_RK2.
9
14
15use mom_checksum_packages, only : mom_thermo_chksum, mom_state_chksum, mom_accel_chksum
16use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
17use mom_cpu_clock, only : clock_component, clock_subcomponent
18use mom_cpu_clock, only : clock_module_driver, clock_module, clock_routine
20use mom_diag_mediator, only : disable_averaging, post_data, safe_alloc_ptr
25use mom_domains, only : to_south, to_west, to_all, cgrid_ne, scalar_pair
26use mom_domains, only : to_north, to_east, omit_corners
27use mom_domains, only : create_group_pass, do_group_pass, group_pass_type
28use mom_domains, only : start_group_pass, complete_group_pass, pass_var, pass_vector
29use mom_debugging, only : hchksum, uvchksum, query_debugging_checks
30use mom_error_handler, only : mom_error, mom_mesg, fatal, warning, is_root_pe
33use mom_file_parser, only : get_param, log_version, param_file_type
34use mom_get_input, only : directories
35use mom_io, only : vardesc, var_desc, east_face, north_face
36use mom_restart, only : register_restart_field, register_restart_pair
37use mom_restart, only : query_initialized, set_initialized, save_restart
38use mom_restart, only : only_read_from_restarts
40use mom_time_manager, only : time_type, real_to_time, operator(+)
41use mom_time_manager, only : operator(-), operator(>), operator(*), operator(/)
42
48use mom_continuity, only : continuity, continuity_cs
49use mom_continuity, only : continuity_init, continuity_stencil
52use mom_debugging, only : check_redundant
53use mom_grid, only : ocean_grid_type
58use mom_interface_heights, only : thickness_to_dz, find_col_avg_spv
60use mom_meke_types, only : meke_type
69use mom_self_attr_load, only : sal_cs
80
81implicit none ; private
82
83#include <MOM_memory.h>
84
85!> MOM_dynamics_split_RK2b module control structure
86type, public :: mom_dyn_split_rk2b_cs ; private
87 real allocable_, dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: &
88 cau, & !< CAu = f*v - u.grad(u) [L T-2 ~> m s-2]
89 cau_pred, & !< The predictor step value of CAu = f*v - u.grad(u) [L T-2 ~> m s-2]
90 pfu, & !< PFu = -dM/dx [L T-2 ~> m s-2]
91 pfu_stokes, & !< PFu_Stokes = -d/dx int_r (u_L*duS/dr) [L T-2 ~> m s-2]
92 diffu !< Zonal acceleration due to convergence of the along-isopycnal stress tensor [L T-2 ~> m s-2]
93
94 real allocable_, dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: &
95 cav, & !< CAv = -f*u - u.grad(v) [L T-2 ~> m s-2]
96 cav_pred, & !< The predictor step value of CAv = -f*u - u.grad(v) [L T-2 ~> m s-2]
97 pfv, & !< PFv = -dM/dy [L T-2 ~> m s-2]
98 pfv_stokes, & !< PFv_Stokes = -d/dy int_r (v_L*dvS/dr) [L T-2 ~> m s-2]
99 diffv !< Meridional acceleration due to convergence of the along-isopycnal stress tensor [L T-2 ~> m s-2]
100
101 real allocable_, dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: visc_rem_u
102 !< Both the fraction of the zonal momentum originally in a
103 !! layer that remains after a time-step of viscosity, and the
104 !! fraction of a time-step worth of a barotropic acceleration
105 !! that a layer experiences after viscosity is applied [nondim].
106 !! Nondimensional between 0 (at the bottom) and 1 (far above).
107 real allocable_, dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: u_accel_bt
108 !< The zonal layer accelerations due to the difference between
109 !! the barotropic accelerations and the baroclinic accelerations
110 !! that were fed into the barotopic calculation [L T-2 ~> m s-2]
111 real allocable_, dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: visc_rem_v
112 !< Both the fraction of the meridional momentum originally in
113 !! a layer that remains after a time-step of viscosity, and the
114 !! fraction of a time-step worth of a barotropic acceleration
115 !! that a layer experiences after viscosity is applied [nondim].
116 !! Nondimensional between 0 (at the bottom) and 1 (far above).
117 real allocable_, dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: v_accel_bt
118 !< The meridional layer accelerations due to the difference between
119 !! the barotropic accelerations and the baroclinic accelerations
120 !! that were fed into the barotopic calculation [L T-2 ~> m s-2]
121
122 ! The following variables are only used with the split time stepping scheme.
123 real allocable_, dimension(NIMEM_,NJMEM_) :: eta !< Instantaneous free surface height (in Boussinesq
124 !! mode) or column mass anomaly (in non-Boussinesq
125 !! mode) [H ~> m or kg m-2]
126 real allocable_, dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: u_av !< layer x-velocity with vertical mean replaced by
127 !! time-mean barotropic velocity over a baroclinic
128 !! timestep [L T-1 ~> m s-1]
129 real allocable_, dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: v_av !< layer y-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_,NJMEM_,NKMEM_) :: h_av !< arithmetic mean of two successive layer
133 !! thicknesses [H ~> m or kg m-2]
134 real allocable_, dimension(NIMEM_,NJMEM_) :: eta_pf !< instantaneous SSH used in calculating PFu and
135 !! PFv [H ~> m or kg m-2]
136 real allocable_, dimension(NIMEMB_PTR_,NJMEM_) :: uhbt !< average x-volume or mass flux determined by the
137 !! barotropic solver [H L2 T-1 ~> m3 s-1 or kg s-1].
138 !! uhbt is roughly equal to the vertical sum of uh.
139 real allocable_, dimension(NIMEM_,NJMEMB_PTR_) :: vhbt !< average y-volume or mass flux determined by the
140 !! barotropic solver [H L2 T-1 ~> m3 s-1 or kg s-1].
141 !! vhbt is roughly equal to vertical sum of vh.
142 real allocable_, dimension(NIMEM_,NJMEM_,NKMEM_) :: pbce !< pbce times eta gives the baroclinic pressure
143 !! anomaly in each layer due to free surface height
144 !! anomalies [L2 H-1 T-2 ~> m s-2 or m4 kg-1 s-2].
145 real allocable_, dimension(NIMEMB_PTR_,NJMEM_) :: du_av_inst !< The barotropic zonal velocity increment
146 !! between filtered and instantaneous velocities
147 !! [L T-1 ~> m s-1]
148 real allocable_, dimension(NIMEM_,NJMEMB_PTR_) :: dv_av_inst !< The barotropic meridional velocity increment
149 !! between filtered and instantaneous velocities
150 !! [L T-1 ~> m s-1]
151
152 real, pointer, dimension(:,:) :: taux_bot => null() !< frictional x-bottom stress from the ocean
153 !! to the seafloor [R L Z T-2 ~> Pa]
154 real, pointer, dimension(:,:) :: tauy_bot => null() !< frictional y-bottom stress from the ocean
155 !! to the seafloor [R L Z T-2 ~> Pa]
156 type(bt_cont_type), pointer :: bt_cont => null() !< A structure with elements that describe the
157 !! effective summed open face areas as a function
158 !! of barotropic flow.
159
160 logical :: bt_adj_corr_mass_src !< If true, recalculates the barotropic mass source after
161 !! predictor step. This should make little difference in the
162 !! deep ocean but appears to help for vanished layers.
163 logical :: split_bottom_stress !< If true, provide the bottom stress
164 !! calculated by the vertical viscosity to the
165 !! barotropic solver.
166 logical :: dtbt_use_bt_cont !< If true, use BT_cont to calculate DTBT.
167 logical :: calculate_sal !< If true, calculate self-attraction and loading.
168 logical :: use_tides !< If true, tidal forcing is enabled.
169 logical :: use_ha !< If true, perform inline harmonic analysis.
170 logical :: remap_aux !< If true, apply ALE remapping to all of the auxiliary 3-D
171 !! variables that are needed to reproduce across restarts,
172 !! similarly to what is done with the primary state variables.
173
174 real :: be !< A nondimensional number from 0.5 to 1 that controls
175 !! the backward weighting of the time stepping scheme [nondim]
176 real :: begw !< A nondimensional number from 0 to 1 that controls
177 !! the extent to which the treatment of gravity waves
178 !! is forward-backward (0) or simulated backward
179 !! Euler (1) [nondim]. 0 is often used.
180 logical :: debug !< If true, write verbose checksums for debugging purposes.
181 logical :: debug_obc !< If true, do additional calls resetting values to help verify the correctness
182 !! of the open boundary condition code.
183 logical :: fpmix = .false. !< If true, applies profiles of momentum flux magnitude and direction.
184 logical :: module_is_initialized = .false. !< Record whether this module has been initialized.
185 logical :: visc_rem_dt_bug = .true. !< If true, recover a bug that uses dt_pred rather than dt for vertvisc_rem
186 !! at the end of predictor.
187
188 !>@{ Diagnostic IDs
189 ! integer :: id_uold = -1, id_vold = -1
190 integer :: id_uh = -1, id_vh = -1
191 integer :: id_umo = -1, id_vmo = -1
192 integer :: id_umo_2d = -1, id_vmo_2d = -1
193 integer :: id_pfu = -1, id_pfv = -1
194 integer :: id_cau = -1, id_cav = -1
195 integer :: id_ueffa = -1, id_veffa = -1
196 ! integer :: id_hf_PFu = -1, id_hf_PFv = -1
197 integer :: id_h_pfu = -1, id_h_pfv = -1
198 integer :: id_hf_pfu_2d = -1, id_hf_pfv_2d = -1
199 integer :: id_intz_pfu_2d = -1, id_intz_pfv_2d = -1
200 integer :: id_pfu_visc_rem = -1, id_pfv_visc_rem = -1
201 ! integer :: id_hf_CAu = -1, id_hf_CAv = -1
202 integer :: id_h_cau = -1, id_h_cav = -1
203 integer :: id_hf_cau_2d = -1, id_hf_cav_2d = -1
204 integer :: id_intz_cau_2d = -1, id_intz_cav_2d = -1
205 integer :: id_cau_visc_rem = -1, id_cav_visc_rem = -1
206 integer :: id_deta_dt = -1
207
208 ! Split scheme only.
209 integer :: id_uav = -1, id_vav = -1
210 integer :: id_u_bt_accel = -1, id_v_bt_accel = -1
211 ! integer :: id_hf_u_BT_accel = -1, id_hf_v_BT_accel = -1
212 integer :: id_h_u_bt_accel = -1, id_h_v_bt_accel = -1
213 integer :: id_hf_u_bt_accel_2d = -1, id_hf_v_bt_accel_2d = -1
214 integer :: id_intz_u_bt_accel_2d = -1, id_intz_v_bt_accel_2d = -1
215 integer :: id_u_bt_accel_visc_rem = -1, id_v_bt_accel_visc_rem = -1
216 !>@}
217
218 type(diag_ctrl), pointer :: diag => null() !< A structure that is used to regulate the
219 !! timing of diagnostic output.
220 type(accel_diag_ptrs), pointer :: adp => null() !< A structure pointing to the various
221 !! accelerations in the momentum equations,
222 !! which can later be used to calculate
223 !! derived diagnostics like energy budgets.
224 type(accel_diag_ptrs), pointer :: ad_pred => null() !< A structure pointing to the various
225 !! predictor step accelerations in the momentum equations,
226 !! which can be used to debug truncations.
227 type(cont_diag_ptrs), pointer :: cdp => null() !< A structure with pointers to various
228 !! terms in the continuity equations,
229 !! which can later be used to calculate
230 !! derived diagnostics like energy budgets.
231
232 ! The remainder of the structure points to child subroutines' control structures.
233 !> A pointer to the horizontal viscosity control structure
234 type(hor_visc_cs) :: hor_visc
235 !> A pointer to the continuity control structure
236 type(continuity_cs) :: continuity_csp
237 !> The CoriolisAdv control structure
238 type(coriolisadv_cs) :: coriolisadv
239 !> A pointer to the PressureForce control structure
240 type(pressureforce_cs) :: pressureforce_csp
241 !> A pointer to a structure containing interface height diffusivities
242 type(vertvisc_cs), pointer :: vertvisc_csp => null()
243 !> A pointer to the set_visc control structure
244 type(set_visc_cs), pointer :: set_visc_csp => null()
245 !> A pointer to the barotropic stepping control structure
246 type(barotropic_cs) :: barotropic_csp
247 !> A pointer to the SAL control structure
248 type(sal_cs) :: sal_csp
249 !> A pointer to the tidal forcing control structure
250 type(tidal_forcing_cs) :: tides_csp
251 !> A pointer to the harmonic analysis control structure
252 type(harmonic_analysis_cs) :: ha_csp
253 !> A pointer to the ALE control structure.
254 type(ale_cs), pointer :: ale_csp => null()
255
256 type(ocean_obc_type), pointer :: obc => null() !< A pointer to an open boundary
257 !! condition type that specifies whether, where, and what open boundary
258 !! conditions are used. If no open BCs are used, this pointer stays
259 !! nullified. Flather OBCs use open boundary_CS as well.
260 !> A pointer to the update_OBC control structure
261 type(update_obc_cs), pointer :: update_obc_csp => null()
262
263 type(group_pass_type) :: pass_eta !< Structure for group halo pass
264 type(group_pass_type) :: pass_visc_rem !< Structure for group halo pass
265 type(group_pass_type) :: pass_uvp !< Structure for group halo pass
266 type(group_pass_type) :: pass_uv_inst !< Structure for group halo pass
267 type(group_pass_type) :: pass_hp_uv !< Structure for group halo pass
268 type(group_pass_type) :: pass_hp_uhvh !< Structure for group halo pass
269 type(group_pass_type) :: pass_h_uv !< Structure for group halo pass
270
272
273
279
280!>@{ CPU time clock IDs
281integer :: id_clock_cor, id_clock_pres, id_clock_vertvisc
282integer :: id_clock_horvisc, id_clock_mom_update
283integer :: id_clock_continuity, id_clock_thick_diff
284integer :: id_clock_btstep, id_clock_btcalc, id_clock_btforce
285integer :: id_clock_pass
286!>@}
287
288contains
289
290!> RK2 splitting for time stepping MOM adiabatic dynamics
291subroutine step_mom_dyn_split_rk2b(u_av, v_av, h, tv, visc, Time_local, dt, forces, &
292 p_surf_begin, p_surf_end, uh, vh, uhtr, vhtr, eta_av, &
293 G, GV, US, CS, calc_dtbt, VarMix, MEKE, thickness_diffuse_CSp, pbv, Waves)
294 type(ocean_grid_type), intent(inout) :: g !< Ocean grid structure
295 type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure
296 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
297 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
298 target, intent(inout) :: u_av !< Zonal velocity [L T-1 ~> m s-1]
299 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
300 target, intent(inout) :: v_av !< Meridional velocity [L T-1 ~> m s-1]
301 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
302 intent(inout) :: h !< Layer thickness [H ~> m or kg m-2]
303 type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic type
304 type(vertvisc_type), intent(inout) :: visc !< Vertical visc, bottom drag, and related
305 type(time_type), intent(in) :: time_local !< Model time at end of time step
306 real, intent(in) :: dt !< Baroclinic dynamics time step [T ~> s]
307 type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces
308 real, dimension(:,:), pointer :: p_surf_begin !< Surface pressure at the start of this dynamic
309 !! time step [R L2 T-2 ~> Pa]
310 real, dimension(:,:), pointer :: p_surf_end !< Surface pressure at the end of this dynamic
311 !! time step [R L2 T-2 ~> Pa]
312 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
313 target, intent(inout) :: uh !< Zonal volume or mass transport
314 !! [H L2 T-1 ~> m3 s-1 or kg s-1]
315 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
316 target, intent(inout) :: vh !< Meridional volume or mass transport
317 !! [H L2 T-1 ~> m3 s-1 or kg s-1]
318 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
319 intent(inout) :: uhtr !< Accumulated zonal volume or mass transport
320 !! since last tracer advection [H L2 ~> m3 or kg]
321 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
322 intent(inout) :: vhtr !< Accumulated meridional volume or mass transport
323 !! since last tracer advection [H L2 ~> m3 or kg]
324 real, dimension(SZI_(G),SZJ_(G)), intent(out) :: eta_av !< Free surface height or column mass
325 !! averaged over time step [H ~> m or kg m-2]
326 type(mom_dyn_split_rk2b_cs), pointer :: cs !< Module control structure
327 logical, intent(in) :: calc_dtbt !< If true, recalculate the barotropic time step
328 type(varmix_cs), intent(inout) :: varmix !< Variable mixing control structure
329 type(meke_type), intent(inout) :: meke !< MEKE fields
330 type(thickness_diffuse_cs), intent(inout) :: thickness_diffuse_csp !< Pointer to a structure containing
331 !! interface height diffusivities
332 type(porous_barrier_type), intent(in) :: pbv !< porous barrier fractional cell metrics
333 type(wave_parameters_cs), optional, pointer :: waves !< A pointer to a structure containing
334 !! fields related to the surface wave conditions
335
336 ! local variables
337 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: up ! Predicted zonal velocity [L T-1 ~> m s-1].
338 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: vp ! Predicted meridional velocity [L T-1 ~> m s-1].
339 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: hp ! Predicted thickness [H ~> m or kg m-2].
340
341 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: dz ! Distance between the interfaces around a layer [Z ~> m]
342 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: ueffa ! Effective Area of U-Faces [H L ~> m2]
343 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: veffa ! Effective Area of V-Faces [H L ~> m2]
344 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: u_bc_accel ! The summed zonal baroclinic accelerations
345 ! of each layer calculated by the non-barotropic
346 ! part of the model [L T-2 ~> m s-2]
347 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: v_bc_accel ! The summed meridional baroclinic accelerations
348 ! of each layer calculated by the non-barotropic
349 ! part of the model [L T-2 ~> m s-2]
350
351 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), target :: u_inst ! Instantaneous zonal velocity [L T-1 ~> m s-1].
352 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), target :: v_inst ! Instantaneous meridional velocity [L T-1 ~> m s-1].
353 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_av ! The average of the layer thicknesses at the beginning
354 ! and end of a time step [H ~> m or kg m-2]
355
356 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), target :: uh_in ! The zonal mass transports that would be
357 ! obtained using the initial velocities [H L2 T-1 ~> m3 s-1 or kg s-1]
358 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), target :: vh_in ! The meridional mass transports that would be
359 ! obtained using the initial velocities [H L2 T-1 ~> m3 s-1 or kg s-1]
360
361 real, dimension(SZI_(G),SZJ_(G)) :: eta_pred ! The predictor value of the free surface height
362 ! or column mass [H ~> m or kg m-2]
363 real, dimension(SZI_(G),SZJ_(G)) :: spv_avg ! The column averaged specific volume [R-1 ~> m3 kg-1]
364 real, dimension(SZI_(G),SZJ_(G)) :: deta_dt ! A diagnostic of the time derivative of the free surface
365 ! height or column mass [H T-1 ~> m s-1 or kg m-2 s-1]
366
367 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: u_old_rad_obc ! The starting zonal velocities, which are
368 ! saved for use in the Flather open boundary condition code [L T-1 ~> m s-1]
369 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: v_old_rad_obc ! The starting meridional velocities, which are
370 ! saved for use in the Flather open boundary condition code [L T-1 ~> m s-1]
371
372 ! GMM, TODO: make these allocatable?
373 ! real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: uold ! u-velocity before vert_visc is applied, for fpmix
374 ! ! [L T-1 ~> m s-1]
375 ! real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: vold ! v-velocity before vert_visc is applied, for fpmix
376 ! ! [L T-1 ~> m s-1]
377 real :: pres_to_eta ! A factor that converts pressures to the units of eta
378 ! [H T2 R-1 L-2 ~> m Pa-1 or kg m-2 Pa-1]
379 real, pointer, dimension(:,:) :: &
380 p_surf => null(), & ! A pointer to the surface pressure [R L2 T-2 ~> Pa]
381 eta_pf_start => null(), & ! The value of eta that corresponds to the starting pressure
382 ! for the barotropic solver [H ~> m or kg m-2]
383 taux_bot => null(), & ! A pointer to the zonal bottom stress in some cases [R L Z T-2 ~> Pa]
384 tauy_bot => null(), & ! A pointer to the meridional bottom stress in some cases [R L Z T-2 ~> Pa]
385 ! This pointer is just used as shorthand for CS%eta.
386 eta => null() ! A pointer to the instantaneous free surface height (in Boussinesq
387 ! mode) or column mass anomaly (in non-Boussinesq mode) [H ~> m or kg m-2]
388
389 real, pointer, dimension(:,:,:) :: &
390 ! These pointers are used to alter which fields are passed to btstep with various options:
391 u_ptr => null(), & ! A pointer to a zonal velocity [L T-1 ~> m s-1]
392 v_ptr => null(), & ! A pointer to a meridional velocity [L T-1 ~> m s-1]
393 uh_ptr => null(), & ! A pointer to a zonal volume or mass transport [H L2 T-1 ~> m3 s-1 or kg s-1]
394 vh_ptr => null() ! A pointer to a meridional volume or mass transport [H L2 T-1 ~> m3 s-1 or kg s-1]
395
396 ! real, dimension(SZI_(G),SZJ_(G)) :: hbl ! Boundary layer depth from Cvmix [H ~> m or kg m-2]
397 real :: dt_pred ! The time step for the predictor part of the baroclinic time stepping [T ~> s].
398 real :: idt_bc ! Inverse of the baroclinic timestep [T-1 ~> s-1]
399 logical :: dyn_p_surf
400 logical :: debug_redundant ! If true, check redundant values on PE boundaries when debugging
401 logical :: bt_cont_bt_thick ! If true, use the BT_cont_type to estimate the
402 ! relative weightings of the layers in calculating
403 ! the barotropic accelerations.
404 !---For group halo pass
405 logical :: showcalltree, sym
406
407 integer :: i, j, k, is, ie, js, je, isq, ieq, jsq, jeq, nz
408 integer :: cont_stencil, obc_stencil
409 integer :: cor_stencil
410
411 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
412 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
413
414 eta => cs%eta
415
416 idt_bc = 1.0 / dt
417
418 sym = g%Domain%symmetric ! switch to include symmetric domain in checksums
419
420 showcalltree = calltree_showquery()
421 if (showcalltree) call calltree_enter("step_MOM_dyn_split_RK2b(), MOM_dynamics_split_RK2b.F90")
422
423 ! Fill in some halo points for arrays that will have halo updates.
424 hp(:,:,:) = h(:,:,:)
425 up(:,:,:) = 0.0 ; vp(:,:,:) = 0.0 ; u_inst(:,:,:) = 0.0 ; v_inst(:,:,:) = 0.0
426
427 ! Update CFL truncation value as function of time
428 call updatecfltruncationvalue(time_local, cs%vertvisc_CSp, us)
429
430 if (cs%debug) then
431 call query_debugging_checks(do_redundant=debug_redundant)
432 call mom_state_chksum("Start predictor ", u_av, v_av, h, uh, vh, g, gv, us, symmetric=sym)
433 if (debug_redundant) then
434 call check_redundant("Start predictor u ", u_av, v_av, g, unscale=us%L_T_to_m_s)
435 call check_redundant("Start predictor uh ", uh, vh, g, unscale=gv%H_to_MKS*us%L_to_m**2*us%s_to_T)
436 endif
437 endif
438
439 dyn_p_surf = associated(p_surf_begin) .and. associated(p_surf_end)
440 if (dyn_p_surf) then
441 p_surf => p_surf_end
442 call safe_alloc_ptr(eta_pf_start,g%isd,g%ied,g%jsd,g%jed)
443 eta_pf_start(:,:) = 0.0
444 else
445 p_surf => forces%p_surf
446 endif
447
448 if (associated(cs%OBC)) then
449 if (cs%debug_OBC) call open_boundary_test_extern_h(g, gv, cs%OBC, h)
450
451 ! Update OBC ramp value as function of time
452 call update_obc_ramp(time_local, cs%OBC, us)
453
454 do k=1,nz ; do j=g%jsd,g%jed ; do i=g%IsdB,g%IedB
455 u_old_rad_obc(i,j,k) = u_av(i,j,k)
456 enddo ; enddo ; enddo
457 do k=1,nz ; do j=g%JsdB,g%JedB ; do i=g%isd,g%ied
458 v_old_rad_obc(i,j,k) = v_av(i,j,k)
459 enddo ; enddo ; enddo
460 endif
461
462 bt_cont_bt_thick = .false.
463 if (associated(cs%BT_cont)) bt_cont_bt_thick = &
464 (allocated(cs%BT_cont%h_u) .and. allocated(cs%BT_cont%h_v))
465
466 if (cs%split_bottom_stress) then
467 taux_bot => cs%taux_bot ; tauy_bot => cs%tauy_bot
468 endif
469
470 !--- begin set up for group halo pass
471
472 cont_stencil = continuity_stencil(cs%continuity_CSp)
473 obc_stencil = 2
474 cor_stencil = coriolisadv_stencil(cs%CoriolisAdv)
475 if (associated(cs%OBC)) then
476 if (cs%OBC%oblique_BCs_exist_globally) obc_stencil = 3
477 endif
478 call cpu_clock_begin(id_clock_pass)
479 call create_group_pass(cs%pass_eta, eta, g%Domain, halo=1)
480 call create_group_pass(cs%pass_visc_rem, cs%visc_rem_u, cs%visc_rem_v, g%Domain, &
481 to_all+scalar_pair, cgrid_ne, halo=max(1,cont_stencil))
482 call create_group_pass(cs%pass_uvp, up, vp, g%Domain, halo=max(1,cont_stencil))
483 call create_group_pass(cs%pass_uv_inst, u_inst, v_inst, g%Domain, halo=max(2,cont_stencil))
484
485 call create_group_pass(cs%pass_hp_uv, hp, g%Domain, halo=cor_stencil)
486 call create_group_pass(cs%pass_hp_uv, u_av, v_av, g%Domain, halo=max(cor_stencil,obc_stencil))
487 call create_group_pass(cs%pass_hp_uv, uh(:,:,:), vh(:,:,:), g%Domain, halo=max(cor_stencil,obc_stencil))
488
489 call create_group_pass(cs%pass_hp_uhvh, hp, g%Domain, halo=cor_stencil)
490 call create_group_pass(cs%pass_hp_uhvh, uh(:,:,:), vh(:,:,:), g%Domain, halo=max(cor_stencil,obc_stencil))
491 if (cor_stencil > 2) then
492 call create_group_pass(cs%pass_hp_uhvh, u_av, v_av, g%Domain, halo=max(cor_stencil,obc_stencil))
493 call create_group_pass(cs%pass_hp_uhvh, h, g%Domain, halo=cor_stencil)
494 endif
495
496 call create_group_pass(cs%pass_h_uv, h, g%Domain, halo=max(cor_stencil,cont_stencil))
497 call create_group_pass(cs%pass_h_uv, u_av, v_av, g%Domain, halo=max(cor_stencil,obc_stencil))
498 call create_group_pass(cs%pass_h_uv, uh(:,:,:), vh(:,:,:), g%Domain, halo=max(cor_stencil,obc_stencil))
499
500 call cpu_clock_end(id_clock_pass)
501 !--- end set up for group halo pass
502
503 ! This calculates the transports and averaged thicknesses that will be used for the
504 ! predictor version of the Coriolis scheme.
505 call cpu_clock_begin(id_clock_continuity)
506 call continuity(u_av, v_av, h, hp, uh, vh, dt, g, gv, us, cs%continuity_CSp, cs%OBC, pbv)
507 call cpu_clock_end(id_clock_continuity)
508
509 if (g%nonblocking_updates) &
510 call start_group_pass(cs%pass_hp_uhvh, g%Domain, clock=id_clock_pass)
511
512! PFu = d/dx M(h,T,S)
513! pbce = dM/deta
514 if (cs%begw == 0.0) call enable_averages(dt, time_local, cs%diag)
515 call cpu_clock_begin(id_clock_pres)
516 call pressureforce(h, tv, cs%PFu, cs%PFv, g, gv, us, cs%PressureForce_CSp, &
517 cs%ALE_CSp, cs%ADp, p_surf, cs%pbce, cs%eta_PF)
518 if (dyn_p_surf) then
519 pres_to_eta = 1.0 / (gv%g_Earth * gv%H_to_RZ)
520 !$OMP parallel do default(shared)
521 do j=jsq,jeq+1 ; do i=isq,ieq+1
522 eta_pf_start(i,j) = cs%eta_PF(i,j) - pres_to_eta * (p_surf_begin(i,j) - p_surf_end(i,j))
523 enddo ; enddo
524 endif
525 ! Stokes shear force contribution to pressure gradient
526 if (present(waves)) then ; if (associated(waves)) then ; if (waves%Stokes_PGF) then
527 call thickness_to_dz(h, tv, dz, g, gv, us, halo_size=1)
528 call stokes_pgf(g, gv, us, dz, u_av, v_av, cs%PFu_Stokes, cs%PFv_Stokes, waves)
529
530 ! We are adding Stokes_PGF to hydrostatic PGF here. The diag PFu/PFv
531 ! will therefore report the sum total PGF and we avoid other
532 ! modifications in the code. The PFu_Stokes is output within the waves routines.
533 if (.not.waves%Passive_Stokes_PGF) then
534 do k=1,nz ; do j=js,je ; do i=isq,ieq
535 cs%PFu(i,j,k) = cs%PFu(i,j,k) + cs%PFu_Stokes(i,j,k)
536 enddo ; enddo ; enddo
537 do k=1,nz ; do j=jsq,jeq ; do i=is,ie
538 cs%PFv(i,j,k) = cs%PFv(i,j,k) + cs%PFv_Stokes(i,j,k)
539 enddo ; enddo ; enddo
540 endif
541 endif ; endif ; endif
542 call cpu_clock_end(id_clock_pres)
543 call disable_averaging(cs%diag)
544 if (showcalltree) call calltree_waypoint("done with PressureForce (step_MOM_dyn_split_RK2b)")
545
546 if (associated(cs%OBC)) then ; if (cs%OBC%update_OBC) then
547 call update_obc_data(cs%OBC, g, gv, us, tv, h, cs%update_OBC_CSp, time_local)
548 endif ; endif
549 if (associated(cs%OBC) .and. cs%debug_OBC) &
550 call open_boundary_zero_normal_flow(cs%OBC, g, gv, cs%PFu, cs%PFv)
551
552 if (g%nonblocking_updates) then
553 call complete_group_pass(cs%pass_hp_uhvh, g%Domain, clock=id_clock_pass)
554 call start_group_pass(cs%pass_eta, g%Domain, clock=id_clock_pass)
555 else
556 call do_group_pass(cs%pass_hp_uhvh, g%Domain, clock=id_clock_pass)
557 endif
558
559 !$OMP parallel do default(shared)
560 do k=1,nz ; do j=js-cor_stencil,je+cor_stencil ; do i=is-cor_stencil,ie+cor_stencil
561 h_av(i,j,k) = 0.5*(h(i,j,k) + hp(i,j,k))
562 enddo ; enddo ; enddo
563
564 ! Calculate a predictor-step estimate of the Coriolis and momentum advection terms
565 ! and horizontal viscous accelerations.
566
567! CAu = -(f+zeta_av)/h_av vh + d/dx KE_av
568 call cpu_clock_begin(id_clock_cor)
569 call coradcalc(u_av, v_av, h_av, uh, vh, cs%CAu_pred, cs%CAv_pred, cs%OBC, cs%AD_pred, &
570 g, gv, us, cs%CoriolisAdv, pbv, waves=waves)
571 call cpu_clock_end(id_clock_cor)
572 if (showcalltree) call calltree_waypoint("done with predictor CorAdCalc (step_MOM_dyn_split_RK2b)")
573
574! diffu = horizontal viscosity terms (u_av)
575 call cpu_clock_begin(id_clock_horvisc)
576 call horizontal_viscosity(u_av, v_av, h_av, uh, vh, cs%diffu, cs%diffv, meke, varmix, g, gv, us, cs%hor_visc, &
577 tv, dt, obc=cs%OBC, bt=cs%barotropic_CSp, td=thickness_diffuse_csp, adp=cs%AD_pred)
578 call cpu_clock_end(id_clock_horvisc)
579 if (showcalltree) call calltree_waypoint("done with predictor horizontal_viscosity (step_MOM_dyn_split_RK2b)")
580
581! u_bc_accel = CAu + PFu + diffu(u[n-1])
582 call cpu_clock_begin(id_clock_btforce)
583 !$OMP parallel do default(shared)
584 do k=1,nz
585 do j=js,je ; do i=isq,ieq
586 u_bc_accel(i,j,k) = (cs%CAu_pred(i,j,k) + cs%PFu(i,j,k)) + cs%diffu(i,j,k)
587 enddo ; enddo
588 do j=jsq,jeq ; do i=is,ie
589 v_bc_accel(i,j,k) = (cs%CAv_pred(i,j,k) + cs%PFv(i,j,k)) + cs%diffv(i,j,k)
590 enddo ; enddo
591 enddo
592 if (associated(cs%OBC)) then
593 call open_boundary_zero_normal_flow(cs%OBC, g, gv, u_bc_accel, v_bc_accel)
594 endif
595 call cpu_clock_end(id_clock_btforce)
596
597 if (cs%debug) then
598 call mom_accel_chksum("pre-btstep accel", cs%CAu_pred, cs%CAv_pred, cs%PFu, cs%PFv, &
599 cs%diffu, cs%diffv, g, gv, us, cs%pbce, u_bc_accel, v_bc_accel, &
600 symmetric=sym)
601 if (debug_redundant) then
602 call check_redundant("pre-btstep CS%CA ", cs%CAu_pred, cs%CAv_pred, g, unscale=us%L_T2_to_m_s2)
603 call check_redundant("pre-btstep CS%PF ", cs%PFu, cs%PFv, g, unscale=us%L_T2_to_m_s2)
604 call check_redundant("pre-btstep CS%diff ", cs%diffu, cs%diffv, g, unscale=us%L_T2_to_m_s2)
605 call check_redundant("pre-btstep u_bc_accel ", u_bc_accel, v_bc_accel, g, unscale=us%L_T2_to_m_s2)
606 endif
607 endif
608
609 call cpu_clock_begin(id_clock_vertvisc)
610 !$OMP parallel do default(shared)
611 do k=1,nz
612 do j=js,je ; do i=isq,ieq
613 up(i,j,k) = g%mask2dCu(i,j) * (u_av(i,j,k) + dt * u_bc_accel(i,j,k))
614 enddo ; enddo
615 do j=jsq,jeq ; do i=is,ie
616 vp(i,j,k) = g%mask2dCv(i,j) * (v_av(i,j,k) + dt * v_bc_accel(i,j,k))
617 enddo ; enddo
618 enddo
619
620 call enable_averages(dt, time_local, cs%diag)
621 call set_viscous_ml(u_av, v_av, h, tv, forces, visc, dt, g, gv, us, cs%set_visc_CSp)
622 call disable_averaging(cs%diag)
623
624 if (cs%debug) then
625 call uvchksum("before vertvisc: up", up, vp, g%HI, haloshift=0, symmetric=sym, unscale=us%L_T_to_m_s)
626 endif
627 call thickness_to_dz(h, tv, dz, g, gv, us, halo_size=1)
628 call vertvisc_coef(up, vp, h, dz, forces, visc, tv, dt, g, gv, us, cs%vertvisc_CSp, cs%OBC, varmix)
629 call vertvisc_remnant(visc, cs%visc_rem_u, cs%visc_rem_v, dt, g, gv, us, cs%vertvisc_CSp)
630 call cpu_clock_end(id_clock_vertvisc)
631 if (showcalltree) call calltree_waypoint("done with vertvisc_coef (step_MOM_dyn_split_RK2b)")
632
633
634 call cpu_clock_begin(id_clock_pass)
635 if (g%nonblocking_updates) then
636 call complete_group_pass(cs%pass_eta, g%Domain)
637 call start_group_pass(cs%pass_visc_rem, g%Domain)
638 else
639 call do_group_pass(cs%pass_eta, g%Domain)
640 call do_group_pass(cs%pass_visc_rem, g%Domain)
641 endif
642 call cpu_clock_end(id_clock_pass)
643
644 call cpu_clock_begin(id_clock_btcalc)
645 ! Calculate the relative layer weights for determining barotropic quantities.
646 if (.not.bt_cont_bt_thick) &
647 call btcalc(h, g, gv, cs%barotropic_CSp, obc=cs%OBC)
648 call bt_mass_source(h, eta, .true., g, gv, cs%barotropic_CSp)
649
650 spv_avg(:,:) = 0.0
651 if ((.not.gv%Boussinesq) .and. associated(cs%OBC)) then
652 ! Determine the column average specific volume if it is needed due to the
653 ! use of Flather open boundary conditions in non-Boussinesq mode.
654 if (open_boundary_query(cs%OBC, apply_flather_obc=.true.)) &
655 call find_col_avg_spv(h, spv_avg, tv, g, gv, us)
656 endif
657 call cpu_clock_end(id_clock_btcalc)
658
659 if (g%nonblocking_updates) &
660 call complete_group_pass(cs%pass_visc_rem, g%Domain, clock=id_clock_pass)
661
662 ! Reconstruct u_inst and v_inst from u_av and v_av.
663 call cpu_clock_begin(id_clock_mom_update)
664 do k=1,nz ; do j=js,je ; do i=isq,ieq
665 u_inst(i,j,k) = u_av(i,j,k) - cs%du_av_inst(i,j) * cs%visc_rem_u(i,j,k)
666 enddo ; enddo ; enddo
667 do k=1,nz ; do j=jsq,jeq ; do i=is,ie
668 v_inst(i,j,k) = v_av(i,j,k) - cs%dv_av_inst(i,j) * cs%visc_rem_v(i,j,k)
669 enddo ; enddo ; enddo
670 call cpu_clock_end(id_clock_mom_update)
671 call do_group_pass(cs%pass_uv_inst, g%Domain, clock=id_clock_pass)
672
673! u_accel_bt = layer accelerations due to barotropic solver
674 call cpu_clock_begin(id_clock_continuity)
675 call continuity(u_inst, v_inst, h, hp, uh_in, vh_in, dt, g, gv, us, cs%continuity_CSp, cs%OBC, pbv, &
676 visc_rem_u=cs%visc_rem_u, visc_rem_v=cs%visc_rem_v, bt_cont=cs%BT_cont)
677 call cpu_clock_end(id_clock_continuity)
678 if (bt_cont_bt_thick) then
679 call btcalc(h, g, gv, cs%barotropic_CSp, cs%BT_cont%h_u, cs%BT_cont%h_v, &
680 obc=cs%OBC)
681 endif
682 if (showcalltree) call calltree_waypoint("done with continuity[BT_cont] (step_MOM_dyn_split_RK2b)")
683
684 uh_ptr => uh_in ; vh_ptr => vh_in ; u_ptr => u_inst ; v_ptr => v_inst
685
686 call cpu_clock_begin(id_clock_btstep)
687 if (calc_dtbt) then
688 if (cs%dtbt_use_bt_cont .and. associated(cs%BT_cont)) then
689 call set_dtbt(g, gv, us, cs%barotropic_CSp, cs%pbce, bt_cont=cs%BT_cont, &
690 time=time_local - real_to_time(dt, unscale=us%T_to_s))
691 else
692 ! In the following call, eta is only used when NONLINEAR_BT_CONTINUITY is True. Otherwise, dtbt is effectively
693 ! calculated with eta=0. Note that NONLINEAR_BT_CONTINUITY is False if BT_CONT is used, which is the default.
694 call set_dtbt(g, gv, us, cs%barotropic_CSp, cs%pbce, eta=eta, &
695 time=time_local - real_to_time(dt, unscale=us%T_to_s))
696 endif
697 endif
698 if (showcalltree) call calltree_enter("btstep(), MOM_barotropic.F90")
699 ! This is the predictor step call to btstep.
700 ! The CS%ADp argument here stores the weights for certain integrated diagnostics.
701 call btstep(u_inst, v_inst, eta, dt, u_bc_accel, v_bc_accel, forces, cs%pbce, cs%eta_PF, u_av, v_av, &
702 cs%u_accel_bt, cs%v_accel_bt, eta_pred, cs%uhbt, cs%vhbt, g, gv, us, &
703 cs%barotropic_CSp, cs%visc_rem_u, cs%visc_rem_v, spv_avg, cs%ADp, cs%OBC, cs%BT_cont, &
704 eta_pf_start, taux_bot, tauy_bot, uh_ptr, vh_ptr, u_ptr, v_ptr)
705 if (showcalltree) call calltree_leave("btstep()")
706 call cpu_clock_end(id_clock_btstep)
707
708 dt_pred = dt * cs%be
709 call cpu_clock_begin(id_clock_mom_update)
710
711! up = u + dt_pred*( u_bc_accel + u_accel_bt )
712 !$OMP parallel do default(shared)
713 do k=1,nz
714 do j=jsq,jeq ; do i=is,ie
715 vp(i,j,k) = g%mask2dCv(i,j) * (v_inst(i,j,k) + dt_pred * &
716 (v_bc_accel(i,j,k) + cs%v_accel_bt(i,j,k)))
717 enddo ; enddo
718 do j=js,je ; do i=isq,ieq
719 up(i,j,k) = g%mask2dCu(i,j) * (u_inst(i,j,k) + dt_pred * &
720 (u_bc_accel(i,j,k) + cs%u_accel_bt(i,j,k)))
721 enddo ; enddo
722 enddo
723 call cpu_clock_end(id_clock_mom_update)
724
725 if (cs%debug) then
726 call uvchksum("Predictor 1 [uv]", up, vp, g%HI, haloshift=0, symmetric=sym, unscale=us%L_T_to_m_s)
727 call hchksum(h, "Predictor 1 h", g%HI, haloshift=1, unscale=gv%H_to_MKS)
728 call uvchksum("Predictor 1 [uv]h", uh, vh, g%HI,haloshift=2, &
729 symmetric=sym, unscale=gv%H_to_MKS*us%L_to_m**2*us%s_to_T)
730! call MOM_state_chksum("Predictor 1", up, vp, h, uh, vh, G, GV, US, haloshift=1)
731 call mom_accel_chksum("Predictor accel", cs%CAu_pred, cs%CAv_pred, cs%PFu, cs%PFv, &
732 cs%diffu, cs%diffv, g, gv, us, cs%pbce, cs%u_accel_bt, cs%v_accel_bt, symmetric=sym)
733 call mom_state_chksum("Predictor 1 init", u_inst, v_inst, h, uh, vh, g, gv, us, haloshift=1, &
734 symmetric=sym)
735 if (debug_redundant) then
736 call check_redundant("Predictor 1 up", up, vp, g, unscale=us%L_T_to_m_s)
737 call check_redundant("Predictor 1 uh", uh, vh, g, unscale=gv%H_to_MKS*us%L_to_m**2*us%s_to_T)
738 endif
739 endif
740
741! up <- up + dt_pred d/dz visc d/dz up
742! u_av <- u_av + dt_pred d/dz visc d/dz u_av
743 call cpu_clock_begin(id_clock_vertvisc)
744 if (cs%debug) then
745 call uvchksum("0 before vertvisc: [uv]p", up, vp, g%HI,haloshift=0, symmetric=sym, unscale=us%L_T_to_m_s)
746 endif
747
748 ! if (CS%fpmix) then
749 ! uold(:,:,:) = 0.0
750 ! vold(:,:,:) = 0.0
751 ! do k=1,nz ; do j=js,je ; do I=Isq,Ieq
752 ! uold(I,j,k) = up(I,j,k)
753 ! enddo ; enddo ; enddo
754 ! do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
755 ! vold(i,J,k) = vp(i,J,k)
756 ! enddo ; enddo ; enddo
757 ! endif
758
759 call thickness_to_dz(h, tv, dz, g, gv, us, halo_size=1)
760 call vertvisc_coef(up, vp, h, dz, forces, visc, tv, dt_pred, g, gv, us, cs%vertvisc_CSp, &
761 cs%OBC, varmix)
762 call vertvisc(up, vp, h, forces, visc, dt_pred, cs%OBC, cs%AD_pred, cs%CDp, g, &
763 gv, us, cs%vertvisc_CSp, cs%taux_bot, cs%tauy_bot, waves=waves)
764
765 ! if (CS%fpmix) then
766 ! hbl(:,:) = 0.0
767 ! if (associated(visc%h_ML)) hbl(:,:) = visc%h_ML(:,:)
768 ! call vertFPmix(up, vp, uold, vold, hbl, h, forces, &
769 ! dt_pred, G, GV, US, CS%vertvisc_CSp, CS%OBC)
770 ! call vertvisc(up, vp, h, forces, visc, dt_pred, CS%OBC, CS%ADp, CS%CDp, G, &
771 ! GV, US, CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot, waves=waves)
772 ! endif
773
774 if (showcalltree) call calltree_waypoint("done with vertvisc (step_MOM_dyn_split_RK2b)")
775 if (g%nonblocking_updates) then
776 call cpu_clock_end(id_clock_vertvisc)
777 call start_group_pass(cs%pass_uvp, g%Domain, clock=id_clock_pass)
778 call cpu_clock_begin(id_clock_vertvisc)
779 endif
780 if (cs%visc_rem_dt_bug) then
781 call vertvisc_remnant(visc, cs%visc_rem_u, cs%visc_rem_v, dt_pred, g, gv, us, cs%vertvisc_CSp)
782 else
783 call vertvisc_remnant(visc, cs%visc_rem_u, cs%visc_rem_v, dt, g, gv, us, cs%vertvisc_CSp)
784 endif
785 call cpu_clock_end(id_clock_vertvisc)
786
787 call do_group_pass(cs%pass_visc_rem, g%Domain, clock=id_clock_pass)
788 if (g%nonblocking_updates) then
789 call complete_group_pass(cs%pass_uvp, g%Domain, clock=id_clock_pass)
790 else
791 call do_group_pass(cs%pass_uvp, g%Domain, clock=id_clock_pass)
792 endif
793
794 ! uh = u_av * h
795 ! hp = h + dt * div . uh
796 call cpu_clock_begin(id_clock_continuity)
797 call continuity(up, vp, h, hp, uh, vh, dt, g, gv, us, cs%continuity_CSp, cs%OBC, pbv, &
798 cs%uhbt, cs%vhbt, cs%visc_rem_u, cs%visc_rem_v, &
799 u_av, v_av, bt_cont=cs%BT_cont)
800 call cpu_clock_end(id_clock_continuity)
801 if (showcalltree) call calltree_waypoint("done with continuity (step_MOM_dyn_split_RK2b)")
802
803 call do_group_pass(cs%pass_hp_uv, g%Domain, clock=id_clock_pass)
804
805 if (associated(cs%OBC)) then
806
807 if (cs%debug) &
808 call uvchksum("Pre OBC avg [uv]", u_av, v_av, g%HI, haloshift=1, symmetric=sym, unscale=us%L_T_to_m_s)
809
810 call radiation_open_bdry_conds(cs%OBC, u_av, u_old_rad_obc, v_av, v_old_rad_obc, g, gv, us, dt_pred)
811
812 if (cs%debug) &
813 call uvchksum("Post OBC avg [uv]", u_av, v_av, g%HI, haloshift=1, symmetric=sym, unscale=us%L_T_to_m_s)
814 endif
815
816 ! h_av = (h + hp)/2
817 !$OMP parallel do default(shared)
818 do k=1,nz ; do j=js-cor_stencil,je+cor_stencil ; do i=is-cor_stencil,ie+cor_stencil
819 h_av(i,j,k) = 0.5*(h(i,j,k) + hp(i,j,k))
820 enddo ; enddo ; enddo
821
822 ! The correction phase of the time step starts here.
823 call enable_averages(dt, time_local, cs%diag)
824
825 ! Calculate a revised estimate of the free-surface height correction to be
826 ! used in the next call to btstep. This call is at this point so that
827 ! hp can be changed if CS%begw /= 0.
828 ! eta_cor = ... (hidden inside CS%barotropic_CSp)
829 if (cs%BT_adj_corr_mass_src) then
830 call cpu_clock_begin(id_clock_btcalc)
831 call bt_mass_source(hp, eta_pred, .false., g, gv, cs%barotropic_CSp)
832 call cpu_clock_end(id_clock_btcalc)
833 endif
834
835 if (cs%begw /= 0.0) then
836 ! hp <- (1-begw)*h_in + begw*hp
837 ! Back up hp to the value it would have had after a time-step of
838 ! begw*dt. hp is not used again until recalculated by continuity.
839 !$OMP parallel do default(shared)
840 do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1
841 hp(i,j,k) = (1.0-cs%begw)*h(i,j,k) + cs%begw*hp(i,j,k)
842 enddo ; enddo ; enddo
843
844 ! PFu = d/dx M(hp,T,S)
845 ! pbce = dM/deta
846 call cpu_clock_begin(id_clock_pres)
847 call pressureforce(hp, tv, cs%PFu, cs%PFv, g, gv, us, cs%PressureForce_CSp, &
848 cs%ALE_CSp, cs%ADp, p_surf, cs%pbce, cs%eta_PF)
849 ! Stokes shear force contribution to pressure gradient
850 if (present(waves)) then ; if (associated(waves)) then ; if (waves%Stokes_PGF) then
851 call thickness_to_dz(h, tv, dz, g, gv, us, halo_size=1)
852 call stokes_pgf(g, gv, us, dz, u_av, v_av, cs%PFu_Stokes, cs%PFv_Stokes, waves)
853 if (.not.waves%Passive_Stokes_PGF) then
854 do k=1,nz ; do j=js,je ; do i=isq,ieq
855 cs%PFu(i,j,k) = cs%PFu(i,j,k) + cs%PFu_Stokes(i,j,k)
856 enddo ; enddo ; enddo
857 do k=1,nz ; do j=jsq,jeq ; do i=is,ie
858 cs%PFv(i,j,k) = cs%PFv(i,j,k) + cs%PFv_Stokes(i,j,k)
859 enddo ; enddo ; enddo
860 endif
861 endif ; endif ; endif
862 call cpu_clock_end(id_clock_pres)
863 if (showcalltree) call calltree_waypoint("done with PressureForce[hp=(1-b).h+b.h] (step_MOM_dyn_split_RK2b)")
864 endif
865
866 if (bt_cont_bt_thick) then
867 call btcalc(h, g, gv, cs%barotropic_CSp, cs%BT_cont%h_u, cs%BT_cont%h_v, &
868 obc=cs%OBC)
869 if (showcalltree) call calltree_waypoint("done with btcalc[BT_cont_BT_thick] (step_MOM_dyn_split_RK2b)")
870 endif
871
872 if (cs%debug) then
873 call mom_state_chksum("Predictor ", up, vp, hp, uh, vh, g, gv, us, symmetric=sym)
874 call uvchksum("Predictor avg [uv]", u_av, v_av, g%HI, haloshift=1, symmetric=sym, unscale=us%L_T_to_m_s)
875 call hchksum(h_av, "Predictor avg h", g%HI, haloshift=2, unscale=gv%H_to_MKS)
876 ! call MOM_state_chksum("Predictor avg ", u_av, v_av, h_av, uh, vh, G, GV, US)
877 if (debug_redundant) then
878 call check_redundant("Predictor up ", up, vp, g, unscale=us%L_T_to_m_s)
879 call check_redundant("Predictor uh ", uh, vh, g, unscale=gv%H_to_MKS*us%L_to_m**2*us%s_to_T)
880 endif
881 endif
882
883! diffu = horizontal viscosity terms (u_av)
884 call cpu_clock_begin(id_clock_horvisc)
885 call horizontal_viscosity(u_av, v_av, h_av, uh, vh, cs%diffu, cs%diffv, meke, varmix, g, gv, us, cs%hor_visc, &
886 tv, dt, obc=cs%OBC, bt=cs%barotropic_CSp, td=thickness_diffuse_csp, adp=cs%ADp)
887 call cpu_clock_end(id_clock_horvisc)
888 if (showcalltree) call calltree_waypoint("done with horizontal_viscosity (step_MOM_dyn_split_RK2b)")
889
890! CAu = -(f+zeta_av)/h_av vh + d/dx KE_av
891 call cpu_clock_begin(id_clock_cor)
892 call coradcalc(u_av, v_av, h_av, uh, vh, cs%CAu, cs%CAv, cs%OBC, cs%ADp, &
893 g, gv, us, cs%CoriolisAdv, pbv, waves=waves)
894 call cpu_clock_end(id_clock_cor)
895 if (showcalltree) call calltree_waypoint("done with CorAdCalc (step_MOM_dyn_split_RK2b)")
896
897! Calculate the momentum forcing terms for the barotropic equations.
898
899! u_bc_accel = CAu + PFu + diffu(u[n-1])
900 call cpu_clock_begin(id_clock_btforce)
901 !$OMP parallel do default(shared)
902 do k=1,nz
903 do j=js,je ; do i=isq,ieq
904 u_bc_accel(i,j,k) = (cs%Cau(i,j,k) + cs%PFu(i,j,k)) + cs%diffu(i,j,k)
905 enddo ; enddo
906 do j=jsq,jeq ; do i=is,ie
907 v_bc_accel(i,j,k) = (cs%Cav(i,j,k) + cs%PFv(i,j,k)) + cs%diffv(i,j,k)
908 enddo ; enddo
909 enddo
910 if (associated(cs%OBC)) then
911 call open_boundary_zero_normal_flow(cs%OBC, g, gv, u_bc_accel, v_bc_accel)
912 endif
913 call cpu_clock_end(id_clock_btforce)
914
915 if (cs%debug) then
916 call mom_accel_chksum("corr pre-btstep accel", cs%CAu, cs%CAv, cs%PFu, cs%PFv, &
917 cs%diffu, cs%diffv, g, gv, us, cs%pbce, u_bc_accel, v_bc_accel, &
918 symmetric=sym)
919 if (debug_redundant) then
920 call check_redundant("corr pre-btstep CS%CA ", cs%CAu, cs%CAv, g, unscale=us%L_T2_to_m_s2)
921 call check_redundant("corr pre-btstep CS%PF ", cs%PFu, cs%PFv, g, unscale=us%L_T2_to_m_s2)
922 call check_redundant("corr pre-btstep CS%diff ", cs%diffu, cs%diffv, g, unscale=us%L_T2_to_m_s2)
923 call check_redundant("corr pre-btstep u_bc_accel ", u_bc_accel, v_bc_accel, g, unscale=us%L_T2_to_m_s2)
924 endif
925 endif
926
927 ! u_accel_bt = layer accelerations due to barotropic solver
928 ! pbce = dM/deta
929 call cpu_clock_begin(id_clock_btstep)
930
931 uh_ptr => uh ; vh_ptr => vh ; u_ptr => u_av ; v_ptr => v_av
932
933 if (showcalltree) call calltree_enter("btstep(), MOM_barotropic.F90")
934 ! This is the corrector step call to btstep.
935 call btstep(u_inst, v_inst, eta, dt, u_bc_accel, v_bc_accel, forces, cs%pbce, cs%eta_PF, u_av, v_av, &
936 cs%u_accel_bt, cs%v_accel_bt, eta_pred, cs%uhbt, cs%vhbt, g, gv, us, &
937 cs%barotropic_CSp, cs%visc_rem_u, cs%visc_rem_v, spv_avg, cs%ADp, cs%OBC, cs%BT_cont, &
938 eta_pf_start, taux_bot, tauy_bot, uh_ptr, vh_ptr, u_ptr, v_ptr, etaav=eta_av)
939 if (cs%id_deta_dt>0) then
940 do j=js,je ; do i=is,ie ; deta_dt(i,j) = (eta_pred(i,j) - eta(i,j))*idt_bc ; enddo ; enddo
941 endif
942 do j=js,je ; do i=is,ie ; eta(i,j) = eta_pred(i,j) ; enddo ; enddo
943
944 call cpu_clock_end(id_clock_btstep)
945 if (showcalltree) call calltree_leave("btstep()")
946
947 if (cs%debug .and. debug_redundant) then
948 call check_redundant("u_accel_bt ", cs%u_accel_bt, cs%v_accel_bt, g, unscale=us%L_T2_to_m_s2)
949 endif
950
951 ! u = u + dt*( u_bc_accel + u_accel_bt )
952 call cpu_clock_begin(id_clock_mom_update)
953 !$OMP parallel do default(shared)
954 do k=1,nz
955 do j=js,je ; do i=isq,ieq
956 u_inst(i,j,k) = g%mask2dCu(i,j) * (u_inst(i,j,k) + dt * &
957 (u_bc_accel(i,j,k) + cs%u_accel_bt(i,j,k)))
958 enddo ; enddo
959 do j=jsq,jeq ; do i=is,ie
960 v_inst(i,j,k) = g%mask2dCv(i,j) * (v_inst(i,j,k) + dt * &
961 (v_bc_accel(i,j,k) + cs%v_accel_bt(i,j,k)))
962 enddo ; enddo
963 enddo
964 call cpu_clock_end(id_clock_mom_update)
965
966 if (cs%debug) then
967 call uvchksum("Corrector 1 [uv]", u_inst, v_inst, g%HI, haloshift=0, symmetric=sym, unscale=us%L_T_to_m_s)
968 call hchksum(h, "Corrector 1 h", g%HI, haloshift=1, unscale=gv%H_to_MKS)
969 call uvchksum("Corrector 1 [uv]h", uh, vh, g%HI, haloshift=2, &
970 symmetric=sym, unscale=gv%H_to_MKS*us%L_to_m**2*us%s_to_T)
971 ! call MOM_state_chksum("Corrector 1", u_inst, v_inst, h, uh, vh, G, GV, US, haloshift=1)
972 call mom_accel_chksum("Corrector accel", cs%CAu, cs%CAv, cs%PFu, cs%PFv, &
973 cs%diffu, cs%diffv, g, gv, us, cs%pbce, cs%u_accel_bt, cs%v_accel_bt, &
974 symmetric=sym)
975 endif
976
977 ! u <- u + dt d/dz visc d/dz u
978 ! u_av <- u_av + dt d/dz visc d/dz u_av
979 call cpu_clock_begin(id_clock_vertvisc)
980
981 ! if (CS%fpmix) then
982 ! uold(:,:,:) = 0.0
983 ! vold(:,:,:) = 0.0
984 ! do k=1,nz ; do j=js,je ; do I=Isq,Ieq
985 ! uold(I,j,k) = u_inst(I,j,k)
986 ! enddo ; enddo ; enddo
987 ! do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
988 ! vold(i,J,k) = v_inst(i,J,k)
989 ! enddo ; enddo ; enddo
990 ! endif
991
992 call thickness_to_dz(h, tv, dz, g, gv, us, halo_size=1)
993 call vertvisc_coef(u_inst, v_inst, h, dz, forces, visc, tv, dt, g, gv, us, cs%vertvisc_CSp, cs%OBC, varmix)
994 call vertvisc(u_inst, v_inst, h, forces, visc, dt, cs%OBC, cs%ADp, cs%CDp, g, gv, us, &
995 cs%vertvisc_CSp, cs%taux_bot, cs%tauy_bot, waves=waves)
996
997 ! if (CS%fpmix) then
998 ! call vertFPmix(u_inst, v_inst, uold, vold, hbl, h, forces, dt, &
999 ! G, GV, US, CS%vertvisc_CSp, CS%OBC)
1000 ! call vertvisc(u_inst, v_inst, h, forces, visc, dt, CS%OBC, CS%ADp, CS%CDp, G, GV, US, &
1001 ! CS%vertvisc_CSp, CS%taux_bot, CS%tauy_bot, waves=waves)
1002 ! endif
1003
1004 if (g%nonblocking_updates) then
1005 call cpu_clock_end(id_clock_vertvisc)
1006 call start_group_pass(cs%pass_uv_inst, g%Domain, clock=id_clock_pass)
1007 call cpu_clock_begin(id_clock_vertvisc)
1008 endif
1009 call vertvisc_remnant(visc, cs%visc_rem_u, cs%visc_rem_v, dt, g, gv, us, cs%vertvisc_CSp)
1010 call cpu_clock_end(id_clock_vertvisc)
1011 if (showcalltree) call calltree_waypoint("done with vertvisc (step_MOM_dyn_split_RK2b)")
1012
1013 call do_group_pass(cs%pass_visc_rem, g%Domain, clock=id_clock_pass)
1014 if (g%nonblocking_updates) then
1015 call complete_group_pass(cs%pass_uv_inst, g%Domain, clock=id_clock_pass)
1016 else
1017 call do_group_pass(cs%pass_uv_inst, g%Domain, clock=id_clock_pass)
1018 endif
1019
1020 ! uh = u_av * h
1021 ! h = h + dt * div . uh
1022 ! u_av and v_av adjusted so their mass transports match uhbt and vhbt.
1023 call cpu_clock_begin(id_clock_continuity)
1024
1025 call continuity(u_inst, v_inst, h, h, uh, vh, dt, g, gv, us, cs%continuity_CSp, cs%OBC, pbv, &
1026 cs%uhbt, cs%vhbt, cs%visc_rem_u, cs%visc_rem_v, u_av, v_av, &
1027 du_cor=cs%du_av_inst, dv_cor=cs%dv_av_inst)
1028
1029 ! This tests the ability to readjust the instantaneous velocity, and here it changes answers only at roundoff.
1030 ! do k=1,nz ; do j=js,je ; do I=Isq,Ieq
1031 ! u_inst(I,j,k) = u_av(I,j,k) - CS%du_av_inst(I,j) * CS%visc_rem_u(I,j,k)
1032 ! enddo ; enddo ; enddo
1033 ! do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
1034 ! v_inst(i,J,k) = v_av(i,J,k) - CS%dv_av_inst(i,J) * CS%visc_rem_v(i,J,k)
1035 ! enddo ; enddo ; enddo
1036
1037 call cpu_clock_end(id_clock_continuity)
1038
1039 call do_group_pass(cs%pass_h_uv, g%Domain, clock=id_clock_pass)
1040
1041 ! Whenever thickness changes let the diag manager know, target grids
1042 ! for vertical remapping may need to be regenerated.
1043 call diag_update_remap_grids(cs%diag)
1044 if (showcalltree) call calltree_waypoint("done with continuity (step_MOM_dyn_split_RK2b)")
1045
1046 if (associated(cs%OBC)) then
1047 call radiation_open_bdry_conds(cs%OBC, u_av, u_old_rad_obc, v_av, v_old_rad_obc, g, gv, us, dt)
1048 endif
1049
1050 !$OMP parallel do default(shared)
1051 do k=1,nz
1052 do j=js-2,je+2 ; do i=isq-2,ieq+2
1053 uhtr(i,j,k) = uhtr(i,j,k) + uh(i,j,k)*dt
1054 enddo ; enddo
1055 do j=jsq-2,jeq+2 ; do i=is-2,ie+2
1056 vhtr(i,j,k) = vhtr(i,j,k) + vh(i,j,k)*dt
1057 enddo ; enddo
1058 enddo
1059
1060 if (associated(cs%OBC)) then
1061 call update_segment_thickness_reservoirs(g, gv, uhtr, vhtr, h, cs%OBC)
1062 endif
1063
1064 ! if (CS%fpmix) then
1065 ! if (CS%id_uold > 0) call post_data(CS%id_uold, uold, CS%diag)
1066 ! if (CS%id_vold > 0) call post_data(CS%id_vold, vold, CS%diag)
1067 ! endif
1068
1069 ! The time-averaged free surface height has already been set by the last call to btstep.
1070
1071 ! Deallocate this memory to avoid a memory leak. ### We should revisit how this array is declared. -RWH
1072 if (dyn_p_surf .and. associated(eta_pf_start)) deallocate(eta_pf_start)
1073
1074 ! Here various terms used in to update the momentum equations are
1075 ! offered for time averaging.
1076 if (cs%id_PFu > 0) call post_data(cs%id_PFu, cs%PFu, cs%diag)
1077 if (cs%id_PFv > 0) call post_data(cs%id_PFv, cs%PFv, cs%diag)
1078 if (cs%id_CAu > 0) call post_data(cs%id_CAu, cs%CAu, cs%diag)
1079 if (cs%id_CAv > 0) call post_data(cs%id_CAv, cs%CAv, cs%diag)
1080
1081 ! Here the thickness fluxes are offered for time averaging.
1082 if (cs%id_uh > 0) call post_data(cs%id_uh, uh, cs%diag)
1083 if (cs%id_vh > 0) call post_data(cs%id_vh, vh, cs%diag)
1084 if (cs%id_uav > 0) call post_data(cs%id_uav, u_av, cs%diag)
1085 if (cs%id_vav > 0) call post_data(cs%id_vav, v_av, cs%diag)
1086 if (cs%id_u_BT_accel > 0) call post_data(cs%id_u_BT_accel, cs%u_accel_bt, cs%diag)
1087 if (cs%id_v_BT_accel > 0) call post_data(cs%id_v_BT_accel, cs%v_accel_bt, cs%diag)
1088
1089 ! Calculate effective areas and post data
1090 if (cs%id_ueffA > 0) then
1091 ueffa(:,:,:) = 0
1092 do k=1,nz ; do j=js,je ; do i=isq,ieq
1093 if (abs(u_av(i,j,k)) > 0.) ueffa(i,j,k) = uh(i,j,k) / u_av(i,j,k)
1094 enddo ; enddo ; enddo
1095 call post_data(cs%id_ueffA, ueffa, cs%diag)
1096 endif
1097
1098 if (cs%id_veffA > 0) then
1099 veffa(:,:,:) = 0
1100 do k=1,nz ; do j=jsq,jeq ; do i=is,ie
1101 if (abs(v_av(i,j,k)) > 0.) veffa(i,j,k) = vh(i,j,k) / v_av(i,j,k)
1102 enddo ; enddo ; enddo
1103 call post_data(cs%id_veffA, veffa, cs%diag)
1104 endif
1105
1106 ! Diagnostics of the fractional thicknesses times momentum budget terms
1107 ! 3D diagnostics hf_PFu etc. are commented because there is no clarity on proper remapping grid option.
1108 ! The code is retained for debugging purposes in the future.
1109 !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)
1110 !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)
1111 !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)
1112 !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)
1113 !if (CS%id_hf_u_BT_accel > 0) &
1114 ! call post_product_u(CS%id_hf_u_BT_accel, CS%u_accel_bt, CS%ADp%diag_hfrac_u, G, nz, CS%diag)
1115 !if (CS%id_hf_v_BT_accel > 0) &
1116 ! call post_product_v(CS%id_hf_v_BT_accel, CS%v_accel_bt, CS%ADp%diag_hfrac_v, G, nz, CS%diag)
1117
1118 ! Diagnostics for the vertical sum of layer thickness x prssure force accelerations
1119 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)
1120 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)
1121
1122 ! Diagnostics for thickness-weighted vertically averaged prssure force accelerations
1123 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)
1124 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)
1125
1126 ! Diagnostics for thickness x prssure force accelerations
1127 if (cs%id_h_PFu > 0) call post_product_u(cs%id_h_PFu, cs%PFu, cs%ADp%diag_hu, g, nz, cs%diag)
1128 if (cs%id_h_PFv > 0) call post_product_v(cs%id_h_PFv, cs%PFv, cs%ADp%diag_hv, g, nz, cs%diag)
1129
1130 ! Diagnostics of Coriolis acceleratations
1131 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)
1132 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)
1133 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)
1134 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)
1135 if (cs%id_h_CAu > 0) call post_product_u(cs%id_h_CAu, cs%CAu, cs%ADp%diag_hu, g, nz, cs%diag)
1136 if (cs%id_h_CAv > 0) call post_product_v(cs%id_h_CAv, cs%CAv, cs%ADp%diag_hv, g, nz, cs%diag)
1137
1138 ! Diagnostics of barotropic solver acceleratations
1139 if (cs%id_intz_u_BT_accel_2d > 0) &
1140 call post_product_sum_u(cs%id_intz_u_BT_accel_2d, cs%u_accel_bt, cs%ADp%diag_hu, g, nz, cs%diag)
1141 if (cs%id_intz_v_BT_accel_2d > 0) &
1142 call post_product_sum_v(cs%id_intz_v_BT_accel_2d, cs%v_accel_bt, cs%ADp%diag_hv, g, nz, cs%diag)
1143 if (cs%id_hf_u_BT_accel_2d > 0) &
1144 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)
1145 if (cs%id_hf_v_BT_accel_2d > 0) &
1146 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)
1147 if (cs%id_h_u_BT_accel > 0) &
1148 call post_product_u(cs%id_h_u_BT_accel, cs%u_accel_bt, cs%ADp%diag_hu, g, nz, cs%diag)
1149 if (cs%id_h_v_BT_accel > 0) &
1150 call post_product_v(cs%id_h_v_BT_accel, cs%v_accel_bt, cs%ADp%diag_hv, g, nz, cs%diag)
1151
1152 ! Diagnostics for momentum budget terms multiplied by visc_rem_[uv],
1153 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)
1154 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)
1155 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)
1156 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)
1157 if (cs%id_u_BT_accel_visc_rem > 0) &
1158 call post_product_u(cs%id_u_BT_accel_visc_rem, cs%u_accel_bt, cs%ADp%visc_rem_u, g, nz, cs%diag)
1159 if (cs%id_v_BT_accel_visc_rem > 0) &
1160 call post_product_v(cs%id_v_BT_accel_visc_rem, cs%v_accel_bt, cs%ADp%visc_rem_v, g, nz, cs%diag)
1161
1162 ! Diagnostics related to changes in eta
1163 if (cs%id_deta_dt > 0) call post_data(cs%id_deta_dt, deta_dt, cs%diag)
1164
1165 if (cs%debug) then
1166 call mom_state_chksum("Corrector ", u_av, v_av, h, uh, vh, g, gv, us, symmetric=sym)
1167 ! call uvchksum("Corrector inst [uv]", u_inst, v_inst, G%HI, symmetric=sym, unscale=US%L_T_to_m_s)
1168 endif
1169
1170 if (showcalltree) call calltree_leave("step_MOM_dyn_split_RK2b()")
1171
1172end subroutine step_mom_dyn_split_rk2b
1173
1174!> This subroutine sets up any auxiliary restart variables that are specific
1175!! to the split-explicit time stepping scheme. All variables registered here should
1176!! have the ability to be recreated if they are not present in a restart file.
1177subroutine register_restarts_dyn_split_rk2b(HI, GV, US, param_file, CS, restart_CS, uh, vh)
1178 type(hor_index_type), intent(in) :: hi !< Horizontal index structure
1179 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
1180 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1181 type(param_file_type), intent(in) :: param_file !< parameter file
1182 type(mom_dyn_split_rk2b_cs), pointer :: cs !< module control structure
1183 type(mom_restart_cs), intent(inout) :: restart_cs !< MOM restart control structure
1184 real, dimension(SZIB_(HI),SZJ_(HI),SZK_(GV)), &
1185 target, intent(inout) :: uh !< zonal volume or mass transport [H L2 T-1 ~> m3 s-1 or kg s-1]
1186 real, dimension(SZI_(HI),SZJB_(HI),SZK_(GV)), &
1187 target, intent(inout) :: vh !< merid volume or mass transport [H L2 T-1 ~> m3 s-1 or kg s-1]
1188
1189 type(vardesc) :: vd(2)
1190 character(len=48) :: thickness_units, flux_units
1191
1192 integer :: isd, ied, jsd, jed, nz, isdb, iedb, jsdb, jedb
1193 character(len=40) :: mdl = "MOM_dynamics_split_RK2b" ! This module's name.
1194
1195 isd = hi%isd ; ied = hi%ied ; jsd = hi%jsd ; jed = hi%jed ; nz = gv%ke
1196 isdb = hi%IsdB ; iedb = hi%IedB ; jsdb = hi%JsdB ; jedb = hi%JedB
1197
1198 ! This is where a control structure specific to this module would be allocated.
1199 if (associated(cs)) then
1200 call mom_error(warning, "register_restarts_dyn_split_RK2b called with an associated "// &
1201 "control structure.")
1202 return
1203 endif
1204 allocate(cs)
1205
1206 alloc_(cs%diffu(isdb:iedb,jsd:jed,nz)) ; cs%diffu(:,:,:) = 0.0
1207 alloc_(cs%diffv(isd:ied,jsdb:jedb,nz)) ; cs%diffv(:,:,:) = 0.0
1208 alloc_(cs%CAu(isdb:iedb,jsd:jed,nz)) ; cs%CAu(:,:,:) = 0.0
1209 alloc_(cs%CAv(isd:ied,jsdb:jedb,nz)) ; cs%CAv(:,:,:) = 0.0
1210 alloc_(cs%CAu_pred(isdb:iedb,jsd:jed,nz)) ; cs%CAu_pred(:,:,:) = 0.0
1211 alloc_(cs%CAv_pred(isd:ied,jsdb:jedb,nz)) ; cs%CAv_pred(:,:,:) = 0.0
1212 alloc_(cs%PFu(isdb:iedb,jsd:jed,nz)) ; cs%PFu(:,:,:) = 0.0
1213 alloc_(cs%PFv(isd:ied,jsdb:jedb,nz)) ; cs%PFv(:,:,:) = 0.0
1214 alloc_(cs%du_av_inst(isdb:iedb,jsd:jed)) ; cs%du_av_inst(:,:) = 0.0
1215 alloc_(cs%dv_av_inst(isd:ied,jsdb:jedb)) ; cs%dv_av_inst(:,:) = 0.0
1216
1217 allocate(cs%taux_bot(isdb:iedb,jsd:jed), source = 0.0)
1218 allocate(cs%tauy_bot(isd:ied,jsdb:jedb), source = 0.0)
1219
1220 alloc_(cs%eta(isd:ied,jsd:jed)) ; cs%eta(:,:) = 0.0
1221
1222 thickness_units = get_thickness_units(gv)
1223 flux_units = get_flux_units(gv)
1224
1225 if (gv%Boussinesq) then
1226 call register_restart_field(cs%eta, "sfc", .false., restart_cs, &
1227 longname="Free surface Height", units=thickness_units, conversion=gv%H_to_mks)
1228 else
1229 call register_restart_field(cs%eta, "p_bot", .false., restart_cs, &
1230 longname="Bottom Pressure", units=thickness_units, conversion=gv%H_to_mks)
1231 endif
1232
1233 ! These are needed to reconstruct the phase in the barotorpic solution.
1234 vd(1) = var_desc("du_avg_inst", "m s-1", &
1235 "Barotropic velocity increment between instantaneous and filtered zonal velocities", 'u', '1')
1236 vd(2) = var_desc("dv_avg_inst", "m s-1", &
1237 "Barotropic velocity increment between instantaneous and filtered meridional velocities", 'v', '1')
1238 call register_restart_pair(cs%du_av_inst, cs%dv_av_inst, vd(1), vd(2), .false., restart_cs, &
1239 conversion=us%L_T_to_m_s)
1240
1241 call register_barotropic_restarts(hi, gv, us, param_file, cs%barotropic_CSp, restart_cs)
1242
1243 call get_param(param_file, mdl, "SPLIT_BOTTOM_STRESS", cs%split_bottom_stress, &
1244 "If true, provide the bottom stress calculated by the "//&
1245 "vertical viscosity to the barotropic solver.", default=.false.,&
1246 do_not_log=.true.)
1247
1248 if (cs%split_bottom_stress) then
1249 vd(1) = var_desc("taux_bot", "kg m-1 s-2", "Zonal bottom stress", 'u', '1')
1250 vd(2) = var_desc("tauy_bot", "kg m-1 s-2", "Meridional bottom stress", 'v', '1')
1251 call register_restart_pair(cs%taux_bot, cs%tauy_bot, vd(1), vd(2), .false., restart_cs, &
1252 conversion=us%RLZ_T2_to_Pa)
1253 endif
1254
1256
1257!> This subroutine does remapping for the auxiliary restart variables that are used
1258!! with the split RK2 time stepping scheme.
1259subroutine remap_dyn_split_rk2b_aux_vars(G, GV, CS, h_old_u, h_old_v, h_new_u, h_new_v, ALE_CSp)
1260 type(ocean_grid_type), intent(inout) :: g !< ocean grid structure
1261 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
1262 type(mom_dyn_split_rk2b_cs), pointer :: cs !< module control structure
1263 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
1264 intent(in) :: h_old_u !< Source grid thickness at zonal
1265 !! velocity points [H ~> m or kg m-2]
1266 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
1267 intent(in) :: h_old_v !< Source grid thickness at meridional
1268 !! velocity points [H ~> m or kg m-2]
1269 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
1270 intent(in) :: h_new_u !< Destination grid thickness at zonal
1271 !! velocity points [H ~> m or kg m-2]
1272 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
1273 intent(in) :: h_new_v !< Destination grid thickness at meridional
1274 !! velocity points [H ~> m or kg m-2]
1275 type(ale_cs), pointer :: ale_csp !< ALE control structure to use when remapping
1276
1277 return
1278
1279end subroutine remap_dyn_split_rk2b_aux_vars
1280
1281!> This subroutine initializes all of the variables that are used by this
1282!! dynamic core, including diagnostics and the cpu clocks.
1283subroutine initialize_dyn_split_rk2b(u, v, h, tv, uh, vh, eta, Time, G, GV, US, param_file, &
1284 diag, CS, HA_CSp, restart_CS, dt, Accel_diag, Cont_diag, MIS, &
1285 VarMix, MEKE, thickness_diffuse_CSp, &
1286 OBC, update_OBC_CSp, ALE_CSp, set_visc, &
1287 visc, dirs, ntrunc, pbv, calc_dtbt, cont_stencil, dyn_h_stencil)
1288 type(ocean_grid_type), intent(inout) :: g !< ocean grid structure
1289 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
1290 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1291 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
1292 intent(inout) :: u !< zonal velocity [L T-1 ~> m s-1]
1293 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
1294 intent(inout) :: v !< merid velocity [L T-1 ~> m s-1]
1295 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1296 intent(inout) :: h !< layer thickness [H ~> m or kg m-2]
1297 type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic type
1298 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
1299 target, intent(inout) :: uh !< zonal volume/mass transport [H L2 T-1 ~> m3 s-1 or kg s-1]
1300 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
1301 target, intent(inout) :: vh !< merid volume/mass transport [H L2 T-1 ~> m3 s-1 or kg s-1]
1302 real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: eta !< free surface height or column mass [H ~> m or kg m-2]
1303 type(time_type), target, intent(in) :: time !< current model time
1304 type(param_file_type), intent(in) :: param_file !< parameter file for parsing
1305 type(diag_ctrl), target, intent(inout) :: diag !< to control diagnostics
1306 type(mom_dyn_split_rk2b_cs), pointer :: cs !< module control structure
1307 type(harmonic_analysis_cs), pointer :: ha_csp !< A pointer to the control structure of the
1308 !! harmonic analysis module
1309 type(mom_restart_cs), intent(inout) :: restart_cs !< MOM restart control structure
1310 real, intent(in) :: dt !< time step [T ~> s]
1311 type(accel_diag_ptrs), target, intent(inout) :: accel_diag !< points to momentum equation terms for
1312 !! budget analysis
1313 type(cont_diag_ptrs), target, intent(inout) :: cont_diag !< points to terms in continuity equation
1314 type(ocean_internal_state), intent(inout) :: mis !< "MOM6 internal state" used to pass
1315 !! diagnostic pointers
1316 type(varmix_cs), intent(inout) :: varmix !< points to spatially variable viscosities
1317 type(meke_type), intent(inout) :: meke !< MEKE fields
1318 type(thickness_diffuse_cs), intent(inout) :: thickness_diffuse_csp !< Pointer to the control structure
1319 !! used for the isopycnal height diffusive transport.
1320 type(ocean_obc_type), pointer :: obc !< points to OBC related fields
1321 type(update_obc_cs), pointer :: update_obc_csp !< points to OBC update related fields
1322 type(ale_cs), pointer :: ale_csp !< points to ALE control structure
1323 type(set_visc_cs), target, intent(in) :: set_visc !< set_visc control structure
1324 type(vertvisc_type), intent(inout) :: visc !< vertical viscosities, bottom drag, and related
1325 type(directories), intent(in) :: dirs !< contains directory paths
1326 integer, target, intent(inout) :: ntrunc !< A target for the variable that records
1327 !! the number of times the velocity is
1328 !! truncated (this should be 0).
1329 logical, intent(out) :: calc_dtbt !< If true, recalculate the barotropic time step
1330 type(porous_barrier_type), intent(in) :: pbv !< porous barrier fractional cell metrics
1331 integer, intent(out) :: cont_stencil !< The stencil for thickness
1332 !! from the continuity solver.
1333 integer, intent(out) :: dyn_h_stencil !< The stencil for thickness for the
1334 !! dynamics based on the continuity
1335 !! solver and Coriolis scheme.
1336
1337 ! local variables
1338 character(len=40) :: mdl = "MOM_dynamics_split_RK2b" ! This module's name.
1339 ! This include declares and sets the variable "version".
1340# include "version_variable.h"
1341 character(len=48) :: thickness_units, flux_units, eta_rest_name
1342 logical :: debug_truncations
1343 logical :: enable_bugs ! If true, the defaults for recently added bug-fix flags are set to
1344 ! recreate the bugs, or if false bugs are only used if actively selected.
1345 logical :: visc_rem_bug ! Stores the value of runtime paramter VISC_REM_BUG.
1346
1347 integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz
1348 integer :: isdb, iedb, jsdb, jedb
1349 integer :: nc ! Number of tidal constituents to be harmonically analyzed
1350 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
1351 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
1352 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
1353
1354 if (.not.associated(cs)) call mom_error(fatal, &
1355 "initialize_dyn_split_RK2b called with an unassociated control structure.")
1356 if (cs%module_is_initialized) then
1357 call mom_error(warning, "initialize_dyn_split_RK2b called with a control "// &
1358 "structure that has already been initialized.")
1359 return
1360 endif
1361 cs%module_is_initialized = .true.
1362
1363 cs%diag => diag
1364
1365 call log_version(param_file, mdl, version, "")
1366 call get_param(param_file, mdl, "TIDES", cs%use_tides, &
1367 "If true, apply tidal momentum forcing.", default=.false.)
1368 call get_param(param_file, mdl, "CALCULATE_SAL", cs%calculate_SAL, &
1369 "If true, calculate self-attraction and loading.", default=cs%use_tides)
1370 call get_param(param_file, mdl, "USE_HA", cs%use_HA, &
1371 "If true, perform inline harmonic analysis.", default=.false.)
1372 call get_param(param_file, mdl, "HA_N_CONST", nc, &
1373 "Number of tidal constituents to be harmonically analyzed.", &
1374 default=0, do_not_log=.not.cs%use_HA)
1375 if (nc<=0) cs%use_HA = .false.
1376 call get_param(param_file, mdl, "BE", cs%be, &
1377 "If SPLIT is true, BE determines the relative weighting "//&
1378 "of a 2nd-order Runga-Kutta baroclinic time stepping "//&
1379 "scheme (0.5) and a backward Euler scheme (1) that is "//&
1380 "used for the Coriolis and inertial terms. BE may be "//&
1381 "from 0.5 to 1, but instability may occur near 0.5. "//&
1382 "BE is also applicable if SPLIT is false and USE_RK2 "//&
1383 "is true.", units="nondim", default=0.6)
1384 call get_param(param_file, mdl, "BEGW", cs%begw, &
1385 "If SPLIT is true, BEGW is a number from 0 to 1 that "//&
1386 "controls the extent to which the treatment of gravity "//&
1387 "waves is forward-backward (0) or simulated backward "//&
1388 "Euler (1). 0 is almost always used. "//&
1389 "If SPLIT is false and USE_RK2 is true, BEGW can be "//&
1390 "between 0 and 0.5 to damp gravity waves.", &
1391 units="nondim", default=0.0)
1392 call get_param(param_file, mdl, "SET_DTBT_USE_BT_CONT", cs%dtbt_use_bt_cont, &
1393 "If true, use BT_CONT to calculate DTBT if possible.", default=.false.)
1394 call get_param(param_file, mdl, "SPLIT_BOTTOM_STRESS", cs%split_bottom_stress, &
1395 "If true, provide the bottom stress calculated by the "//&
1396 "vertical viscosity to the barotropic solver.", default=.false.)
1397 call get_param(param_file, mdl, "BT_ADJ_CORR_MASS_SRC", cs%BT_adj_corr_mass_src, &
1398 "If true, recalculates the barotropic mass source after "//&
1399 "predictor step. This should make little difference in the "//&
1400 "deep ocean but appears to help for vanished layers. If false, "//&
1401 "uses the same mass source as from the predictor step.", default=.true.)
1402 ! call get_param(param_file, mdl, "FPMIX", CS%fpmix, &
1403 ! "If true, apply profiles of momentum flux magnitude and direction.", &
1404 ! default=.false.)
1405 cs%fpmix = .false.
1406 call get_param(param_file, mdl, "REMAP_AUXILIARY_VARS", cs%remap_aux, &
1407 "If true, apply ALE remapping to all of the auxiliary 3-dimensional "//&
1408 "variables that are needed to reproduce across restarts, similarly to "//&
1409 "what is already being done with the primary state variables. "//&
1410 "The default should be changed to true.", default=.false., do_not_log=.true.)
1411 call get_param(param_file, mdl, "DEBUG", cs%debug, &
1412 "If true, write out verbose debugging data.", &
1413 default=.false., debuggingparam=.true.)
1414 call get_param(param_file, mdl, "OBC_DEBUGGING_TESTS", cs%debug_OBC, &
1415 "If true, do additional calls resetting certain values to help verify the "//&
1416 "correctness of the open boundary condition code.", &
1417 default=.false., old_name="DEBUG_OBC", debuggingparam=.true., do_not_log=.true.)
1418 call get_param(param_file, mdl, "DEBUG_TRUNCATIONS", debug_truncations, &
1419 default=.false.)
1420 call get_param(param_file, mdl, "ENABLE_BUGS_BY_DEFAULT", enable_bugs, &
1421 default=.true., do_not_log=.true.) ! This is logged from MOM.F90.
1422 call get_param(param_file, mdl, "VISC_REM_BUG", visc_rem_bug, &
1423 "If true, visc_rem_[uv] in split mode is incorrectly calculated or accounted "//&
1424 "for in two places. This parameter controls the defaults of two individual "//&
1425 "flags, VISC_REM_TIMESTEP_BUG in MOM_dynamics_split_RK2(b) and "//&
1426 "VISC_REM_BT_WEIGHT_BUG in MOM_barotropic.", default=.false.)
1427 call get_param(param_file, mdl, "VISC_REM_TIMESTEP_BUG", cs%visc_rem_dt_bug, &
1428 "If true, recover a bug that uses dt_pred rather than dt in "//&
1429 "vertvisc_remnant() at the end of predictor stage for the following "//&
1430 "continuity() and btstep() calls in the corrector step. Default of this flag "//&
1431 "is set by VISC_REM_BUG", default=visc_rem_bug)
1432
1433
1434 alloc_(cs%uhbt(isdb:iedb,jsd:jed)) ; cs%uhbt(:,:) = 0.0
1435 alloc_(cs%vhbt(isd:ied,jsdb:jedb)) ; cs%vhbt(:,:) = 0.0
1436 alloc_(cs%visc_rem_u(isdb:iedb,jsd:jed,nz)) ; cs%visc_rem_u(:,:,:) = 0.0
1437 alloc_(cs%visc_rem_v(isd:ied,jsdb:jedb,nz)) ; cs%visc_rem_v(:,:,:) = 0.0
1438 alloc_(cs%eta_PF(isd:ied,jsd:jed)) ; cs%eta_PF(:,:) = 0.0
1439 alloc_(cs%pbce(isd:ied,jsd:jed,nz)) ; cs%pbce(:,:,:) = 0.0
1440
1441 alloc_(cs%u_accel_bt(isdb:iedb,jsd:jed,nz)) ; cs%u_accel_bt(:,:,:) = 0.0
1442 alloc_(cs%v_accel_bt(isd:ied,jsdb:jedb,nz)) ; cs%v_accel_bt(:,:,:) = 0.0
1443 alloc_(cs%PFu_Stokes(isdb:iedb,jsd:jed,nz)) ; cs%PFu_Stokes(:,:,:) = 0.0
1444 alloc_(cs%PFv_Stokes(isd:ied,jsdb:jedb,nz)) ; cs%PFv_Stokes(:,:,:) = 0.0
1445
1446 mis%diffu => cs%diffu
1447 mis%diffv => cs%diffv
1448 mis%PFu => cs%PFu
1449 mis%PFv => cs%PFv
1450 mis%CAu => cs%CAu
1451 mis%CAv => cs%CAv
1452 mis%pbce => cs%pbce
1453 mis%u_accel_bt => cs%u_accel_bt
1454 mis%v_accel_bt => cs%v_accel_bt
1455
1456 cs%ADp => accel_diag
1457 cs%CDp => cont_diag
1458 accel_diag%diffu => cs%diffu
1459 accel_diag%diffv => cs%diffv
1460 accel_diag%PFu => cs%PFu
1461 accel_diag%PFv => cs%PFv
1462 accel_diag%CAu => cs%CAu
1463 accel_diag%CAv => cs%CAv
1464 accel_diag%u_accel_bt => cs%u_accel_bt
1465 accel_diag%v_accel_bt => cs%v_accel_bt
1466
1467 allocate(cs%AD_pred)
1468 cs%AD_pred%diffu => cs%diffu
1469 cs%AD_pred%diffv => cs%diffv
1470 cs%AD_pred%PFu => cs%PFu
1471 cs%AD_pred%PFv => cs%PFv
1472 cs%AD_pred%CAu => cs%CAu_pred
1473 cs%AD_pred%CAv => cs%CAv_pred
1474 cs%AD_pred%u_accel_bt => cs%u_accel_bt
1475 cs%AD_pred%v_accel_bt => cs%v_accel_bt
1476
1477! Accel_diag%pbce => CS%pbce
1478! Accel_diag%u_accel_bt => CS%u_accel_bt ; Accel_diag%v_accel_bt => CS%v_accel_bt
1479! Accel_diag%u_av => CS%u_av ; Accel_diag%v_av => CS%v_av
1480
1481 call continuity_init(time, g, gv, us, param_file, diag, cs%continuity_CSp, cs%OBC)
1482 cont_stencil = continuity_stencil(cs%continuity_CSp)
1483 call coriolisadv_init(time, g, gv, us, param_file, diag, cs%ADp, cs%CoriolisAdv)
1484 dyn_h_stencil = max(cont_stencil, coriolisadv_stencil(cs%CoriolisAdv))
1485 if (cs%calculate_SAL) call sal_init(h, tv, g, gv, us, param_file, cs%SAL_CSp, restart_cs)
1486 if (cs%use_tides) call tidal_forcing_init(time, g, us, param_file, cs%tides_CSp)
1487 if (cs%use_HA) then
1488 call ha_init(time, us, param_file, nc, cs%HA_CSp)
1489 ha_csp => cs%HA_CSp
1490 else
1491 ha_csp => null()
1492 endif
1493 call pressureforce_init(time, g, gv, us, param_file, diag, cs%PressureForce_CSp, cs%ADp, &
1494 cs%SAL_CSp, cs%tides_CSp)
1495 call hor_visc_init(time, g, gv, us, param_file, diag, cs%hor_visc, adp=cs%ADp)
1496 call vertvisc_init(mis, time, g, gv, us, param_file, diag, cs%ADp, dirs, &
1497 ntrunc, cs%vertvisc_CSp)
1498 cs%set_visc_CSp => set_visc
1499 call updatecfltruncationvalue(time, cs%vertvisc_CSp, us, &
1500 activate=is_new_run(restart_cs) )
1501
1502 if (associated(ale_csp)) cs%ALE_CSp => ale_csp
1503 if (associated(obc)) then
1504 cs%OBC => obc
1505 if (obc%ramp) call update_obc_ramp(time, cs%OBC, us, &
1506 activate=is_new_run(restart_cs) )
1507 endif
1508 if (associated(update_obc_csp)) cs%update_OBC_CSp => update_obc_csp
1509
1510 eta_rest_name = "sfc" ; if (.not.gv%Boussinesq) eta_rest_name = "p_bot"
1511 if (.not. query_initialized(cs%eta, trim(eta_rest_name), restart_cs)) then
1512 ! Estimate eta based on the layer thicknesses - h. With the Boussinesq
1513 ! approximation, eta is the free surface height anomaly, while without it
1514 ! eta is the mass of ocean per unit area. eta always has the same
1515 ! dimensions as h, either m or kg m-3.
1516 ! CS%eta(:,:) = 0.0 already from initialization.
1517 if (gv%Boussinesq) then
1518 do j=js,je ; do i=is,ie ; cs%eta(i,j) = -gv%Z_to_H * g%bathyT(i,j) ; enddo ; enddo
1519 endif
1520 do k=1,nz ; do j=js,je ; do i=is,ie
1521 cs%eta(i,j) = cs%eta(i,j) + h(i,j,k)
1522 enddo ; enddo ; enddo
1523 call set_initialized(cs%eta, trim(eta_rest_name), restart_cs)
1524 endif
1525 ! Copy eta into an output array.
1526 do j=js,je ; do i=is,ie ; eta(i,j) = cs%eta(i,j) ; enddo ; enddo
1527
1528 call barotropic_init(u, v, h, time, g, gv, us, param_file, diag, &
1529 cs%barotropic_CSp, restart_cs, calc_dtbt, cs%BT_cont, &
1530 cs%OBC, cs%SAL_CSp, ha_csp)
1531
1532 flux_units = get_flux_units(gv)
1533 thickness_units = get_thickness_units(gv)
1534 cs%id_uh = register_diag_field('ocean_model', 'uh', diag%axesCuL, time, &
1535 'Zonal Thickness Flux', flux_units, conversion=gv%H_to_MKS*us%L_to_m**2*us%s_to_T, &
1536 y_cell_method='sum', v_extensive=.true.)
1537 cs%id_vh = register_diag_field('ocean_model', 'vh', diag%axesCvL, time, &
1538 'Meridional Thickness Flux', flux_units, conversion=gv%H_to_MKS*us%L_to_m**2*us%s_to_T, &
1539 x_cell_method='sum', v_extensive=.true.)
1540
1541 cs%id_CAu = register_diag_field('ocean_model', 'CAu', diag%axesCuL, time, &
1542 'Zonal Coriolis and Advective Acceleration', 'm s-2', conversion=us%L_T2_to_m_s2)
1543 cs%id_CAv = register_diag_field('ocean_model', 'CAv', diag%axesCvL, time, &
1544 'Meridional Coriolis and Advective Acceleration', 'm s-2', conversion=us%L_T2_to_m_s2)
1545 cs%id_PFu = register_diag_field('ocean_model', 'PFu', diag%axesCuL, time, &
1546 'Zonal Pressure Force Acceleration', 'm s-2', conversion=us%L_T2_to_m_s2)
1547 cs%id_PFv = register_diag_field('ocean_model', 'PFv', diag%axesCvL, time, &
1548 'Meridional Pressure Force Acceleration', 'm s-2', conversion=us%L_T2_to_m_s2)
1549 cs%id_ueffA = register_diag_field('ocean_model', 'ueffA', diag%axesCuL, time, &
1550 'Effective U-Face Area', 'm^2', conversion=gv%H_to_m*us%L_to_m, &
1551 y_cell_method='sum', v_extensive=.true.)
1552 cs%id_veffA = register_diag_field('ocean_model', 'veffA', diag%axesCvL, time, &
1553 'Effective V-Face Area', 'm^2', conversion=gv%H_to_m*us%L_to_m, &
1554 x_cell_method='sum', v_extensive=.true.)
1555 if (gv%Boussinesq) then
1556 cs%id_deta_dt = register_diag_field('ocean_model', 'deta_dt', diag%axesT1, time, &
1557 'Barotropic SSH tendency due to dynamics', trim(thickness_units)//' s-1', conversion=gv%H_to_MKS*us%s_to_T)
1558 else
1559 cs%id_deta_dt = register_diag_field('ocean_model', 'deta_dt', diag%axesT1, time, &
1560 'Barotropic column-mass tendency due to dynamics', trim(thickness_units)//' s-1', &
1561 conversion=gv%H_to_mks*us%s_to_T)
1562 endif
1563
1564 !CS%id_hf_PFu = register_diag_field('ocean_model', 'hf_PFu', diag%axesCuL, Time, &
1565 ! 'Fractional Thickness-weighted Zonal Pressure Force Acceleration', &
1566 ! 'm s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)
1567 !if (CS%id_hf_PFu > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_u,IsdB,IedB,jsd,jed,nz)
1568
1569 !CS%id_hf_PFv = register_diag_field('ocean_model', 'hf_PFv', diag%axesCvL, Time, &
1570 ! 'Fractional Thickness-weighted Meridional Pressure Force Acceleration', &
1571 ! 'm s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)
1572 !if (CS%id_hf_PFv > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz)
1573
1574 !CS%id_hf_CAu = register_diag_field('ocean_model', 'hf_CAu', diag%axesCuL, Time, &
1575 ! 'Fractional Thickness-weighted Zonal Coriolis and Advective Acceleration', &
1576 ! 'm s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)
1577 !if (CS%id_hf_CAu > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_u,IsdB,IedB,jsd,jed,nz)
1578
1579 !CS%id_hf_CAv = register_diag_field('ocean_model', 'hf_CAv', diag%axesCvL, Time, &
1580 ! 'Fractional Thickness-weighted Meridional Coriolis and Advective Acceleration', &
1581 ! 'm s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)
1582 !if (CS%id_hf_CAv > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz)
1583
1584 cs%id_hf_PFu_2d = register_diag_field('ocean_model', 'hf_PFu_2d', diag%axesCu1, time, &
1585 'Depth-sum Fractional Thickness-weighted Zonal Pressure Force Acceleration', &
1586 'm s-2', conversion=us%L_T2_to_m_s2)
1587 if (cs%id_hf_PFu_2d > 0) call safe_alloc_ptr(cs%ADp%diag_hfrac_u,isdb,iedb,jsd,jed,nz)
1588
1589 cs%id_hf_PFv_2d = register_diag_field('ocean_model', 'hf_PFv_2d', diag%axesCv1, time, &
1590 'Depth-sum Fractional Thickness-weighted Meridional Pressure Force Acceleration', &
1591 'm s-2', conversion=us%L_T2_to_m_s2)
1592 if (cs%id_hf_PFv_2d > 0) call safe_alloc_ptr(cs%ADp%diag_hfrac_v,isd,ied,jsdb,jedb,nz)
1593
1594 cs%id_h_PFu = register_diag_field('ocean_model', 'h_PFu', diag%axesCuL, time, &
1595 'Thickness Multiplied Zonal Pressure Force Acceleration', &
1596 'm2 s-2', conversion=gv%H_to_m*us%L_T2_to_m_s2)
1597 if (cs%id_h_PFu > 0) call safe_alloc_ptr(cs%ADp%diag_hu,isdb,iedb,jsd,jed,nz)
1598
1599 cs%id_h_PFv = register_diag_field('ocean_model', 'h_PFv', diag%axesCvL, time, &
1600 'Thickness Multiplied Meridional Pressure Force Acceleration', &
1601 'm2 s-2', conversion=gv%H_to_m*us%L_T2_to_m_s2)
1602 if (cs%id_h_PFv > 0) call safe_alloc_ptr(cs%ADp%diag_hv,isd,ied,jsdb,jedb,nz)
1603
1604 cs%id_intz_PFu_2d = register_diag_field('ocean_model', 'intz_PFu_2d', diag%axesCu1, time, &
1605 'Depth-integral of Zonal Pressure Force Acceleration', &
1606 'm2 s-2', conversion=gv%H_to_m*us%L_T2_to_m_s2)
1607 if (cs%id_intz_PFu_2d > 0) call safe_alloc_ptr(cs%ADp%diag_hu,isdb,iedb,jsd,jed,nz)
1608
1609 cs%id_intz_PFv_2d = register_diag_field('ocean_model', 'intz_PFv_2d', diag%axesCv1, time, &
1610 'Depth-integral of Meridional Pressure Force Acceleration', &
1611 'm2 s-2', conversion=gv%H_to_m*us%L_T2_to_m_s2)
1612 if (cs%id_intz_PFv_2d > 0) call safe_alloc_ptr(cs%ADp%diag_hv,isd,ied,jsdb,jedb,nz)
1613
1614 cs%id_hf_CAu_2d = register_diag_field('ocean_model', 'hf_CAu_2d', diag%axesCu1, time, &
1615 'Depth-sum Fractional Thickness-weighted Zonal Coriolis and Advective Acceleration', &
1616 'm s-2', conversion=us%L_T2_to_m_s2)
1617 if (cs%id_hf_CAu_2d > 0) call safe_alloc_ptr(cs%ADp%diag_hfrac_u,isdb,iedb,jsd,jed,nz)
1618
1619 cs%id_hf_CAv_2d = register_diag_field('ocean_model', 'hf_CAv_2d', diag%axesCv1, time, &
1620 'Depth-sum Fractional Thickness-weighted Meridional Coriolis and Advective Acceleration', &
1621 'm s-2', conversion=us%L_T2_to_m_s2)
1622 if (cs%id_hf_CAv_2d > 0) call safe_alloc_ptr(cs%ADp%diag_hfrac_v,isd,ied,jsdb,jedb,nz)
1623
1624 cs%id_h_CAu = register_diag_field('ocean_model', 'h_CAu', diag%axesCuL, time, &
1625 'Thickness Multiplied Zonal Coriolis and Advective Acceleration', &
1626 'm2 s-2', conversion=gv%H_to_m*us%L_T2_to_m_s2)
1627 if (cs%id_h_CAu > 0) call safe_alloc_ptr(cs%ADp%diag_hu,isdb,iedb,jsd,jed,nz)
1628
1629 cs%id_h_CAv = register_diag_field('ocean_model', 'h_CAv', diag%axesCvL, time, &
1630 'Thickness Multiplied Meridional Coriolis and Advective Acceleration', &
1631 'm2 s-2', conversion=gv%H_to_m*us%L_T2_to_m_s2)
1632 if (cs%id_h_CAv > 0) call safe_alloc_ptr(cs%ADp%diag_hv,isd,ied,jsdb,jedb,nz)
1633
1634 cs%id_intz_CAu_2d = register_diag_field('ocean_model', 'intz_CAu_2d', diag%axesCu1, time, &
1635 'Depth-integral of Zonal Coriolis and Advective Acceleration', &
1636 'm2 s-2', conversion=gv%H_to_m*us%L_T2_to_m_s2)
1637 if (cs%id_intz_CAu_2d > 0) call safe_alloc_ptr(cs%ADp%diag_hu,isdb,iedb,jsd,jed,nz)
1638
1639 cs%id_intz_CAv_2d = register_diag_field('ocean_model', 'intz_CAv_2d', diag%axesCv1, time, &
1640 'Depth-integral of Meridional Coriolis and Advective Acceleration', &
1641 'm2 s-2', conversion=gv%H_to_m*us%L_T2_to_m_s2)
1642 if (cs%id_intz_CAv_2d > 0) call safe_alloc_ptr(cs%ADp%diag_hv,isd,ied,jsdb,jedb,nz)
1643
1644 cs%id_uav = register_diag_field('ocean_model', 'uav', diag%axesCuL, time, &
1645 'Barotropic-step Averaged Zonal Velocity', 'm s-1', conversion=us%L_T_to_m_s)
1646 cs%id_vav = register_diag_field('ocean_model', 'vav', diag%axesCvL, time, &
1647 'Barotropic-step Averaged Meridional Velocity', 'm s-1', conversion=us%L_T_to_m_s)
1648
1649 cs%id_u_BT_accel = register_diag_field('ocean_model', 'u_BT_accel', diag%axesCuL, time, &
1650 'Barotropic Anomaly Zonal Acceleration', 'm s-2', conversion=us%L_T2_to_m_s2)
1651 cs%id_v_BT_accel = register_diag_field('ocean_model', 'v_BT_accel', diag%axesCvL, time, &
1652 'Barotropic Anomaly Meridional Acceleration', 'm s-2', conversion=us%L_T2_to_m_s2)
1653
1654 !CS%id_hf_u_BT_accel = register_diag_field('ocean_model', 'hf_u_BT_accel', diag%axesCuL, Time, &
1655 ! 'Fractional Thickness-weighted Barotropic Anomaly Zonal Acceleration', &
1656 ! 'm s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)
1657 !if (CS%id_hf_u_BT_accel > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_u,IsdB,IedB,jsd,jed,nz)
1658
1659 !CS%id_hf_v_BT_accel = register_diag_field('ocean_model', 'hf_v_BT_accel', diag%axesCvL, Time, &
1660 ! 'Fractional Thickness-weighted Barotropic Anomaly Meridional Acceleration', &
1661 ! 'm s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)
1662 !if (CS%id_hf_v_BT_accel > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz)
1663
1664 cs%id_hf_u_BT_accel_2d = register_diag_field('ocean_model', 'hf_u_BT_accel_2d', diag%axesCu1, time, &
1665 'Depth-sum Fractional Thickness-weighted Barotropic Anomaly Zonal Acceleration', &
1666 'm s-2', conversion=us%L_T2_to_m_s2)
1667 if (cs%id_hf_u_BT_accel_2d > 0) call safe_alloc_ptr(cs%ADp%diag_hfrac_u,isdb,iedb,jsd,jed,nz)
1668
1669 cs%id_hf_v_BT_accel_2d = register_diag_field('ocean_model', 'hf_v_BT_accel_2d', diag%axesCv1, time, &
1670 'Depth-sum Fractional Thickness-weighted Barotropic Anomaly Meridional Acceleration', &
1671 'm s-2', conversion=us%L_T2_to_m_s2)
1672 if (cs%id_hf_v_BT_accel_2d > 0) call safe_alloc_ptr(cs%ADp%diag_hfrac_v,isd,ied,jsdb,jedb,nz)
1673
1674 cs%id_h_u_BT_accel = register_diag_field('ocean_model', 'h_u_BT_accel', diag%axesCuL, time, &
1675 'Thickness Multiplied Barotropic Anomaly Zonal Acceleration', &
1676 'm2 s-2', conversion=gv%H_to_m*us%L_T2_to_m_s2)
1677 if (cs%id_h_u_BT_accel > 0) call safe_alloc_ptr(cs%ADp%diag_hu,isdb,iedb,jsd,jed,nz)
1678
1679 cs%id_h_v_BT_accel = register_diag_field('ocean_model', 'h_v_BT_accel', diag%axesCvL, time, &
1680 'Thickness Multiplied Barotropic Anomaly Meridional Acceleration', &
1681 'm2 s-2', conversion=gv%H_to_m*us%L_T2_to_m_s2)
1682 if (cs%id_h_v_BT_accel > 0) call safe_alloc_ptr(cs%ADp%diag_hv,isd,ied,jsdb,jedb,nz)
1683
1684 cs%id_intz_u_BT_accel_2d = register_diag_field('ocean_model', 'intz_u_BT_accel_2d', diag%axesCu1, time, &
1685 'Depth-integral of Barotropic Anomaly Zonal Acceleration', &
1686 'm2 s-2', conversion=gv%H_to_m*us%L_T2_to_m_s2)
1687 if (cs%id_intz_u_BT_accel_2d > 0) call safe_alloc_ptr(cs%ADp%diag_hu,isdb,iedb,jsd,jed,nz)
1688
1689 cs%id_intz_v_BT_accel_2d = register_diag_field('ocean_model', 'intz_v_BT_accel_2d', diag%axesCv1, time, &
1690 'Depth-integral of Barotropic Anomaly Meridional Acceleration', &
1691 'm2 s-2', conversion=gv%H_to_m*us%L_T2_to_m_s2)
1692 if (cs%id_intz_v_BT_accel_2d > 0) call safe_alloc_ptr(cs%ADp%diag_hv,isd,ied,jsdb,jedb,nz)
1693
1694 cs%id_PFu_visc_rem = register_diag_field('ocean_model', 'PFu_visc_rem', diag%axesCuL, time, &
1695 'Zonal Pressure Force Acceleration multiplied by the viscous remnant', &
1696 'm s-2', conversion=us%L_T2_to_m_s2)
1697 if (cs%id_PFu_visc_rem > 0) call safe_alloc_ptr(cs%ADp%visc_rem_u,isdb,iedb,jsd,jed,nz)
1698 cs%id_PFv_visc_rem = register_diag_field('ocean_model', 'PFv_visc_rem', diag%axesCvL, time, &
1699 'Meridional Pressure Force Acceleration multiplied by the viscous remnant', &
1700 'm s-2', conversion=us%L_T2_to_m_s2)
1701 if (cs%id_PFv_visc_rem > 0) call safe_alloc_ptr(cs%ADp%visc_rem_v,isd,ied,jsdb,jedb,nz)
1702
1703 cs%id_CAu_visc_rem = register_diag_field('ocean_model', 'CAu_visc_rem', diag%axesCuL, time, &
1704 'Zonal Coriolis and Advective Acceleration multiplied by the viscous remnant', &
1705 'm s-2', conversion=us%L_T2_to_m_s2)
1706 if (cs%id_CAu_visc_rem > 0) call safe_alloc_ptr(cs%ADp%visc_rem_u,isdb,iedb,jsd,jed,nz)
1707 cs%id_CAv_visc_rem = register_diag_field('ocean_model', 'CAv_visc_rem', diag%axesCvL, time, &
1708 'Meridional Coriolis and Advective Acceleration multiplied by the viscous remnant', &
1709 'm s-2', conversion=us%L_T2_to_m_s2)
1710 if (cs%id_CAv_visc_rem > 0) call safe_alloc_ptr(cs%ADp%visc_rem_v,isd,ied,jsdb,jedb,nz)
1711
1712 cs%id_u_BT_accel_visc_rem = register_diag_field('ocean_model', 'u_BT_accel_visc_rem', diag%axesCuL, time, &
1713 'Barotropic Anomaly Zonal Acceleration multiplied by the viscous remnant', &
1714 'm s-2', conversion=us%L_T2_to_m_s2)
1715 if (cs%id_u_BT_accel_visc_rem > 0) call safe_alloc_ptr(cs%ADp%visc_rem_u,isdb,iedb,jsd,jed,nz)
1716 cs%id_v_BT_accel_visc_rem = register_diag_field('ocean_model', 'v_BT_accel_visc_rem', diag%axesCvL, time, &
1717 'Barotropic Anomaly Meridional Acceleration multiplied by the viscous remnant', &
1718 'm s-2', conversion=us%L_T2_to_m_s2)
1719 if (cs%id_v_BT_accel_visc_rem > 0) call safe_alloc_ptr(cs%ADp%visc_rem_v,isd,ied,jsdb,jedb,nz)
1720
1721 id_clock_cor = cpu_clock_id('(Ocean Coriolis & mom advection)', grain=clock_module)
1722 id_clock_continuity = cpu_clock_id('(Ocean continuity equation)', grain=clock_module)
1723 id_clock_pres = cpu_clock_id('(Ocean pressure force)', grain=clock_module)
1724 id_clock_vertvisc = cpu_clock_id('(Ocean vertical viscosity)', grain=clock_module)
1725 id_clock_horvisc = cpu_clock_id('(Ocean horizontal viscosity)', grain=clock_module)
1726 id_clock_mom_update = cpu_clock_id('(Ocean momentum increments)', grain=clock_module)
1727 id_clock_pass = cpu_clock_id('(Ocean message passing)', grain=clock_module)
1728 id_clock_btcalc = cpu_clock_id('(Ocean barotropic mode calc)', grain=clock_module)
1729 id_clock_btstep = cpu_clock_id('(Ocean barotropic mode stepping)', grain=clock_module)
1730 id_clock_btforce = cpu_clock_id('(Ocean barotropic forcing calc)', grain=clock_module)
1731
1732end subroutine initialize_dyn_split_rk2b
1733
1734
1735!> Close the dyn_split_RK2b module
1736subroutine end_dyn_split_rk2b(CS)
1737 type(mom_dyn_split_rk2b_cs), pointer :: cs !< module control structure
1738
1739 call barotropic_end(cs%barotropic_CSp)
1740
1741 call vertvisc_end(cs%vertvisc_CSp)
1742 deallocate(cs%vertvisc_CSp)
1743
1744 call hor_visc_end(cs%hor_visc)
1745 if (cs%calculate_SAL) call sal_end(cs%SAL_CSp)
1746 if (cs%use_tides) call tidal_forcing_end(cs%tides_CSp)
1747 call coriolisadv_end(cs%CoriolisAdv)
1748
1749 dealloc_(cs%diffu) ; dealloc_(cs%diffv)
1750 dealloc_(cs%CAu) ; dealloc_(cs%CAv)
1751 dealloc_(cs%CAu_pred) ; dealloc_(cs%CAv_pred)
1752 dealloc_(cs%PFu) ; dealloc_(cs%PFv)
1753
1754 if (associated(cs%taux_bot)) deallocate(cs%taux_bot)
1755 if (associated(cs%tauy_bot)) deallocate(cs%tauy_bot)
1756 dealloc_(cs%uhbt) ; dealloc_(cs%vhbt)
1757 dealloc_(cs%du_av_inst) ; dealloc_(cs%dv_av_inst)
1758 dealloc_(cs%u_accel_bt) ; dealloc_(cs%v_accel_bt)
1759 dealloc_(cs%visc_rem_u) ; dealloc_(cs%visc_rem_v)
1760
1761 dealloc_(cs%eta) ; dealloc_(cs%eta_PF) ; dealloc_(cs%pbce)
1762
1763 call dealloc_bt_cont_type(cs%BT_cont)
1764 deallocate(cs%AD_pred)
1765
1766 deallocate(cs)
1767end subroutine end_dyn_split_rk2b
1768
1769
1770!> \namespace mom_dynamics_split_rk2b
1771!!
1772!! This file time steps the adiabatic dynamic core by splitting
1773!! between baroclinic and barotropic modes. It uses a pseudo-second order
1774!! Runge-Kutta time stepping scheme for the baroclinic momentum
1775!! equation and a forward-backward coupling between the baroclinic
1776!! momentum and continuity equations. This split time-stepping
1777!! scheme is described in detail in Hallberg (JCP, 1997). Additional
1778!! issues related to exact tracer conservation and how to
1779!! ensure consistency between the barotropic and layered estimates
1780!! of the free surface height are described in Hallberg and
1781!! Adcroft (Ocean Modelling, 2009). This was the time stepping code
1782!! that is used for most GOLD applications, including GFDL's ESM2G
1783!! Earth system model, and all of the examples provided with the
1784!! MOM code (although several of these solutions are routinely
1785!! verified by comparison with the slower unsplit schemes).
1786!!
1787!! The subroutine step_MOM_dyn_split_RK2b actually does the time
1788!! stepping, while register_restarts_dyn_split_RK2b sets the fields
1789!! that are found in a full restart file with this scheme, and
1790!! initialize_dyn_split_RK2b initializes the cpu clocks that are
1791!! used in this module. For largely historical reasons, this module
1792!! does not have its own control structure, but shares the same
1793!! control structure with MOM.F90 and the other MOM_dynamics_...
1794!! modules.
1795
1796end module mom_dynamics_split_rk2b