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