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