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