MOM_dynamics_split_RK2.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> Time step the adiabatic dynamic core of MOM using RK2 method.
7
12
13use mom_checksum_packages, only : mom_thermo_chksum, mom_state_chksum, mom_accel_chksum
14use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
15use mom_cpu_clock, only : clock_component, clock_subcomponent
16use mom_cpu_clock, only : clock_module_driver, clock_module, clock_routine
19use mom_diag_mediator, only : disable_averaging, post_data, safe_alloc_ptr
24use mom_domains, only : to_south, to_west, to_all, cgrid_ne, scalar_pair
25use mom_domains, only : to_north, to_east, omit_corners
26use mom_domains, only : create_group_pass, do_group_pass, group_pass_type
27use mom_domains, only : start_group_pass, complete_group_pass, pass_var, pass_vector
28use mom_debugging, only : hchksum, uvchksum, query_debugging_checks
29use mom_error_handler, only : mom_error, mom_mesg, fatal, warning, is_root_pe
32use mom_file_parser, only : get_param, log_version, param_file_type
33use mom_get_input, only : directories
34use mom_io, only : vardesc, var_desc, east_face, north_face
35use mom_restart, only : register_restart_field, register_restart_pair
36use mom_restart, only : query_initialized, set_initialized, save_restart
37use mom_restart, only : only_read_from_restarts
39use mom_time_manager, only : time_type, real_to_time, operator(+)
40use mom_time_manager, only : operator(-), operator(>), operator(*), operator(/)
41
47use mom_continuity, only : continuity, continuity_cs
48use mom_continuity, only : continuity_init, continuity_stencil
52use mom_debugging, only : check_redundant
54use mom_grid, only : ocean_grid_type
59use mom_interface_heights, only : thickness_to_dz, find_col_avg_spv
61use mom_meke_types, only : meke_type
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 time=time_local - real_to_time(dt, unscale=us%T_to_s))
679 else
680 ! In the following call, eta is only used when NONLINEAR_BT_CONTINUITY is True. Otherwise, dtbt is effectively
681 ! calculated with eta=0. Note that NONLINEAR_BT_CONTINUITY is False if BT_CONT is used, which is the default.
682 call set_dtbt(g, gv, us, cs%barotropic_CSp, cs%pbce, eta=eta, &
683 time=time_local - real_to_time(dt, unscale=us%T_to_s))
684 endif
685 endif
686 if (showcalltree) call calltree_enter("btstep(), MOM_barotropic.F90")
687 ! This is the predictor step call to btstep.
688 ! The CS%ADp argument here stores the weights for certain integrated diagnostics.
689 call btstep(u_inst, v_inst, eta, dt, u_bc_accel, v_bc_accel, forces, cs%pbce, cs%eta_PF, u_av, v_av, &
690 cs%u_accel_bt, cs%v_accel_bt, eta_pred, cs%uhbt, cs%vhbt, g, gv, us, &
691 cs%barotropic_CSp, cs%visc_rem_u, cs%visc_rem_v, spv_avg, cs%ADp, cs%OBC, cs%BT_cont, &
692 eta_pf_start, taux_bot, tauy_bot, uh_ptr, vh_ptr, u_ptr, v_ptr)
693 if (showcalltree) call calltree_leave("btstep()")
694 call cpu_clock_end(id_clock_btstep)
695
696! up = u + dt_pred*( u_bc_accel + u_accel_bt )
697 dt_pred = dt * cs%be
698 call cpu_clock_begin(id_clock_mom_update)
699
700 !$OMP parallel do default(shared)
701 do k=1,nz
702 do j=jsq,jeq ; do i=is,ie
703 vp(i,j,k) = g%mask2dCv(i,j) * (v_inst(i,j,k) + dt_pred * &
704 (v_bc_accel(i,j,k) + cs%v_accel_bt(i,j,k)))
705 enddo ; enddo
706 do j=js,je ; do i=isq,ieq
707 up(i,j,k) = g%mask2dCu(i,j) * (u_inst(i,j,k) + dt_pred * &
708 (u_bc_accel(i,j,k) + cs%u_accel_bt(i,j,k)))
709 enddo ; enddo
710 enddo
711 call cpu_clock_end(id_clock_mom_update)
712
713 if (cs%debug) then
714 call mom_accel_chksum("Predictor accel", cs%CAu_pred, cs%CAv_pred, cs%PFu, cs%PFv, &
715 cs%diffu, cs%diffv, g, gv, us, cs%pbce, cs%u_accel_bt, cs%v_accel_bt, symmetric=sym)
716 call uvchksum("Predictor 1 [uv]", up, vp, g%HI, haloshift=0, symmetric=sym, unscale=us%L_T_to_m_s)
717 call hchksum(h, "Predictor 1 h", g%HI, haloshift=1, unscale=gv%H_to_MKS)
718 call uvchksum("Predictor 1 [uv]h", uh, vh, g%HI,haloshift=2, &
719 symmetric=sym, unscale=gv%H_to_MKS*us%L_to_m**2*us%s_to_T)
720 ! call MOM_state_chksum("Predictor 1", up, vp, h, uh, vh, G, GV, US, haloshift=1)
721 call mom_state_chksum("Predictor 1 init", u_inst, v_inst, h, uh, vh, g, gv, us, haloshift=1, &
722 symmetric=sym)
723 if (debug_redundant) then
724 call check_redundant("Predictor 1 up", up, vp, g, unscale=us%L_T_to_m_s)
725 call check_redundant("Predictor 1 uh", uh, vh, g, unscale=gv%H_to_MKS*us%L_to_m**2*us%s_to_T)
726 endif
727 endif
728
729! up <- up + dt_pred d/dz visc d/dz up
730! u_av <- u_av + dt_pred d/dz visc d/dz u_av
731 call cpu_clock_begin(id_clock_vertvisc)
732 if (cs%debug) then
733 call uvchksum("0 before vertvisc: [uv]p", up, vp, g%HI,haloshift=0, symmetric=sym, unscale=us%L_T_to_m_s)
734 endif
735
736 if (cs%fpmix) then
737 uold(:,:,:) = 0.0
738 vold(:,:,:) = 0.0
739 do k = 1, nz
740 do j = js , je
741 do i = isq, ieq
742 uold(i,j,k) = up(i,j,k)
743 enddo
744 enddo
745 do j = jsq, jeq
746 do i = is, ie
747 vold(i,j,k) = vp(i,j,k)
748 enddo
749 enddo
750 enddo
751 endif
752
753 call thickness_to_dz(h, tv, dz, g, gv, us, halo_size=1)
754 call vertvisc_coef(up, vp, h, dz, forces, visc, tv, dt_pred, g, gv, us, cs%vertvisc_CSp, &
755 cs%OBC, varmix)
756
757 if (cs%fpmix) then
758 hbl(:,:) = 0.0
759 if (ASSOCIATED(cs%KPP_CSp)) call kpp_get_bld(cs%KPP_CSp, hbl, g, us, m_to_bld_units=gv%m_to_H)
760 if (ASSOCIATED(cs%energetic_PBL_CSp)) &
761 call energetic_pbl_get_mld(cs%energetic_PBL_CSp, hbl, g, us, m_to_mld_units=gv%m_to_H)
762
763 ! lFPpost must be false in the predictor step to avoid averaging into the diagnostics
764 lfppost = .false.
765 call vertfpmix(up, vp, uold, vold, hbl, h, forces, dt_pred, lfppost, cs%Cemp_NL, &
766 g, gv, us, cs%vertvisc_CSp, cs%OBC, waves=waves)
767 call vertvisc(up, vp, h, forces, visc, dt_pred, cs%OBC, cs%AD_pred, cs%CDp, g, &
768 gv, us, cs%vertvisc_CSp, cs%taux_bot, cs%tauy_bot, fpmix=cs%fpmix, waves=waves)
769 else
770 call vertvisc(up, vp, h, forces, visc, dt_pred, cs%OBC, cs%AD_pred, cs%CDp, g, &
771 gv, us, cs%vertvisc_CSp, cs%taux_bot, cs%tauy_bot, waves=waves)
772 endif
773
774 if (showcalltree) call calltree_waypoint("done with vertvisc (step_MOM_dyn_split_RK2)")
775 if (g%nonblocking_updates) then
776 call cpu_clock_end(id_clock_vertvisc)
777 call start_group_pass(cs%pass_uvp, g%Domain, clock=id_clock_pass)
778 call cpu_clock_begin(id_clock_vertvisc)
779 endif
780 if (cs%visc_rem_dt_bug) then
781 call vertvisc_remnant(visc, cs%visc_rem_u, cs%visc_rem_v, dt_pred, g, gv, us, cs%vertvisc_CSp)
782 else
783 call vertvisc_remnant(visc, cs%visc_rem_u, cs%visc_rem_v, dt, g, gv, us, cs%vertvisc_CSp)
784 endif
785 call cpu_clock_end(id_clock_vertvisc)
786
787 call do_group_pass(cs%pass_visc_rem, g%Domain, clock=id_clock_pass)
788 if (g%nonblocking_updates) then
789 call complete_group_pass(cs%pass_uvp, g%Domain, clock=id_clock_pass)
790 else
791 call do_group_pass(cs%pass_uvp, g%Domain, clock=id_clock_pass)
792 endif
793
794 ! uh = u_av * h
795 ! hp = h + dt * div . uh
796 call cpu_clock_begin(id_clock_continuity)
797 call continuity(up, vp, h, hp, uh, vh, dt, g, gv, us, cs%continuity_CSp, cs%OBC, pbv, &
798 uhbt=cs%uhbt, vhbt=cs%vhbt, visc_rem_u=cs%visc_rem_u, visc_rem_v=cs%visc_rem_v, &
799 u_cor=u_av, v_cor=v_av, bt_cont=cs%BT_cont)
800 call cpu_clock_end(id_clock_continuity)
801 if (showcalltree) call calltree_waypoint("done with continuity (step_MOM_dyn_split_RK2)")
802
803 call do_group_pass(cs%pass_hp_uv, g%Domain, clock=id_clock_pass)
804
805 if (associated(cs%OBC)) then
806
807 if (cs%debug) &
808 call uvchksum("Pre OBC avg [uv]", u_av, v_av, g%HI, haloshift=1, symmetric=sym, unscale=us%L_T_to_m_s)
809
810 call radiation_open_bdry_conds(cs%OBC, u_av, u_old_rad_obc, v_av, v_old_rad_obc, g, gv, us, dt_pred)
811
812 if (cs%debug) &
813 call uvchksum("Post OBC avg [uv]", u_av, v_av, g%HI, haloshift=1, symmetric=sym, unscale=us%L_T_to_m_s)
814
815 ! These should be done with a pass that excludes uh & vh.
816! call do_group_pass(CS%pass_hp_uv, G%Domain, clock=id_clock_pass)
817 call pass_vector(u_av, v_av, g%Domain, halo=max(cor_stencil,vel_stencil), clock=id_clock_pass)
818 endif
819
820 ! h_av = (h + hp)/2
821 !$OMP parallel do default(shared)
822 do k=1,nz ; do j=js-cor_stencil,je+cor_stencil ; do i=is-cor_stencil,ie+cor_stencil
823 h_av(i,j,k) = 0.5*(h(i,j,k) + hp(i,j,k))
824 enddo ; enddo ; enddo
825
826 ! The correction phase of the time step starts here.
827 call enable_averages(dt, time_local, cs%diag)
828
829 ! Calculate a revised estimate of the free-surface height correction to be
830 ! used in the next call to btstep. This call is at this point so that
831 ! hp can be changed if CS%begw /= 0.
832 ! eta_cor = ... (hidden inside CS%barotropic_CSp)
833 if (cs%BT_adj_corr_mass_src) then
834 call cpu_clock_begin(id_clock_btcalc)
835 call bt_mass_source(hp, eta_pred, .false., g, gv, cs%barotropic_CSp)
836 call cpu_clock_end(id_clock_btcalc)
837 endif
838
839 if (cs%begw /= 0.0) then
840 ! hp <- (1-begw)*h_in + begw*hp
841 ! Back up hp to the value it would have had after a time-step of
842 ! begw*dt. hp is not used again until recalculated by continuity.
843 !$OMP parallel do default(shared)
844 do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1
845 hp(i,j,k) = (1.0-cs%begw)*h(i,j,k) + cs%begw*hp(i,j,k)
846 enddo ; enddo ; enddo
847
848 ! PFu = d/dx M(hp,T,S)
849 ! pbce = dM/deta
850 call cpu_clock_begin(id_clock_pres)
851 call pressureforce(hp, tv, cs%PFu, cs%PFv, g, gv, us, cs%PressureForce_CSp, &
852 cs%ALE_CSp, cs%ADp, p_surf, cs%pbce, cs%eta_PF)
853 ! Stokes shear force contribution to pressure gradient
854 use_stokes_pgf = present(waves)
855 if (use_stokes_pgf) then
856 use_stokes_pgf = associated(waves)
857 if (use_stokes_pgf) use_stokes_pgf = waves%Stokes_PGF
858 if (use_stokes_pgf) then
859 call thickness_to_dz(h, tv, dz, g, gv, us, halo_size=1)
860 call stokes_pgf(g, gv, us, dz, u_inst, v_inst, cs%PFu_Stokes, cs%PFv_Stokes, waves)
861 if (.not.waves%Passive_Stokes_PGF) then
862 do k=1,nz
863 do j=js,je ; do i=isq,ieq
864 cs%PFu(i,j,k) = cs%PFu(i,j,k) + cs%PFu_Stokes(i,j,k)
865 enddo ; enddo
866 enddo
867 do k=1,nz
868 do j=jsq,jeq ; do i=is,ie
869 cs%PFv(i,j,k) = cs%PFv(i,j,k) + cs%PFv_Stokes(i,j,k)
870 enddo ; enddo
871 enddo
872 endif
873 endif
874 endif
875 call cpu_clock_end(id_clock_pres)
876 if (showcalltree) call calltree_waypoint("done with PressureForce[hp=(1-b).h+b.h] (step_MOM_dyn_split_RK2)")
877 endif
878
879 if (bt_cont_bt_thick) then
880 call btcalc(h, g, gv, cs%barotropic_CSp, cs%BT_cont%h_u, cs%BT_cont%h_v, &
881 obc=cs%OBC)
882 if (showcalltree) call calltree_waypoint("done with btcalc[BT_cont_BT_thick] (step_MOM_dyn_split_RK2)")
883 endif
884
885 if (cs%debug) then
886 call mom_state_chksum("Predictor ", up, vp, hp, uh, vh, g, gv, us, symmetric=sym)
887 call uvchksum("Predictor avg [uv]", u_av, v_av, g%HI, haloshift=1, symmetric=sym, unscale=us%L_T_to_m_s)
888 call hchksum(h_av, "Predictor avg h", g%HI, haloshift=2, unscale=gv%H_to_MKS)
889 ! call MOM_state_chksum("Predictor avg ", u_av, v_av, h_av, uh, vh, G, GV, US)
890 if (debug_redundant) then
891 call check_redundant("Predictor up ", up, vp, g, unscale=us%L_T_to_m_s)
892 call check_redundant("Predictor uh ", uh, vh, g, unscale=gv%H_to_MKS*us%L_to_m**2*us%s_to_T)
893 endif
894 endif
895
896! diffu = horizontal viscosity terms (u_av)
897 call cpu_clock_begin(id_clock_horvisc)
898 call horizontal_viscosity(u_av, v_av, h_av, uh, vh, cs%diffu, cs%diffv, &
899 meke, varmix, g, gv, us, cs%hor_visc, tv, dt, &
900 obc=cs%OBC, bt=cs%barotropic_CSp, td=thickness_diffuse_csp, &
901 adp=cs%ADp, hu_cont=cs%BT_cont%h_u, hv_cont=cs%BT_cont%h_v, stoch=stoch)
902 call cpu_clock_end(id_clock_horvisc)
903 if (showcalltree) call calltree_waypoint("done with horizontal_viscosity (step_MOM_dyn_split_RK2)")
904
905! CAu = -(f+zeta_av)/h_av vh + d/dx KE_av
906 call cpu_clock_begin(id_clock_cor)
907 call coradcalc(u_av, v_av, h_av, uh, vh, cs%CAu, cs%CAv, cs%OBC, cs%ADp, &
908 g, gv, us, cs%CoriolisAdv, pbv, waves=waves)
909 call cpu_clock_end(id_clock_cor)
910 if (showcalltree) call calltree_waypoint("done with CorAdCalc (step_MOM_dyn_split_RK2)")
911
912! Calculate the momentum forcing terms for the barotropic equations.
913
914! u_bc_accel = CAu + PFu + diffu(u[n-1])
915 call cpu_clock_begin(id_clock_btforce)
916 !$OMP parallel do default(shared)
917 do k=1,nz
918 do j=js,je ; do i=isq,ieq
919 u_bc_accel(i,j,k) = (cs%Cau(i,j,k) + cs%PFu(i,j,k)) + cs%diffu(i,j,k)
920 enddo ; enddo
921 do j=jsq,jeq ; do i=is,ie
922 v_bc_accel(i,j,k) = (cs%Cav(i,j,k) + cs%PFv(i,j,k)) + cs%diffv(i,j,k)
923 enddo ; enddo
924 enddo
925 if (associated(cs%OBC)) then
926 call open_boundary_zero_normal_flow(cs%OBC, g, gv, u_bc_accel, v_bc_accel)
927 endif
928 call cpu_clock_end(id_clock_btforce)
929
930 if (cs%debug) then
931 call mom_accel_chksum("corr pre-btstep accel", cs%CAu, cs%CAv, cs%PFu, cs%PFv, &
932 cs%diffu, cs%diffv, g, gv, us, cs%pbce, u_bc_accel, v_bc_accel, &
933 symmetric=sym)
934 if (debug_redundant) then
935 call check_redundant("corr pre-btstep CS%CA ", cs%CAu, cs%CAv, g, unscale=us%L_T2_to_m_s2)
936 call check_redundant("corr pre-btstep CS%PF ", cs%PFu, cs%PFv, g, unscale=us%L_T2_to_m_s2)
937 call check_redundant("corr pre-btstep CS%diff ", cs%diffu, cs%diffv, g, unscale=us%L_T2_to_m_s2)
938 call check_redundant("corr pre-btstep u_bc_accel ", u_bc_accel, v_bc_accel, g, unscale=us%L_T2_to_m_s2)
939 endif
940 endif
941
942 ! u_accel_bt = layer accelerations due to barotropic solver
943 ! pbce = dM/deta
944 call cpu_clock_begin(id_clock_btstep)
945 if (cs%BT_use_layer_fluxes) then
946 uh_ptr => uh ; vh_ptr => vh ; u_ptr => u_av ; v_ptr => v_av
947 endif
948
949 if (showcalltree) call calltree_enter("btstep(), MOM_barotropic.F90")
950 ! This is the corrector step call to btstep.
951 call btstep(u_inst, v_inst, eta, dt, u_bc_accel, v_bc_accel, forces, cs%pbce, cs%eta_PF, u_av, v_av, &
952 cs%u_accel_bt, cs%v_accel_bt, eta_pred, cs%uhbt, cs%vhbt, g, gv, us, &
953 cs%barotropic_CSp, cs%visc_rem_u, cs%visc_rem_v, spv_avg, cs%ADp, cs%OBC, cs%BT_cont, &
954 eta_pf_start, taux_bot, tauy_bot, uh_ptr, vh_ptr, u_ptr, v_ptr, etaav=eta_av)
955 if (cs%id_deta_dt>0) then
956 do j=js,je ; do i=is,ie ; deta_dt(i,j) = (eta_pred(i,j) - eta(i,j))*idt_bc ; enddo ; enddo
957 endif
958 do j=js,je ; do i=is,ie ; eta(i,j) = eta_pred(i,j) ; enddo ; enddo
959
960 call cpu_clock_end(id_clock_btstep)
961 if (showcalltree) call calltree_leave("btstep()")
962
963 if (cs%debug .and. debug_redundant) then
964 call check_redundant("u_accel_bt ", cs%u_accel_bt, cs%v_accel_bt, g, unscale=us%L_T2_to_m_s2)
965 endif
966
967 ! u = u + dt*( u_bc_accel + u_accel_bt )
968 call cpu_clock_begin(id_clock_mom_update)
969 !$OMP parallel do default(shared)
970 do k=1,nz
971 do j=js,je ; do i=isq,ieq
972 u_inst(i,j,k) = g%mask2dCu(i,j) * (u_inst(i,j,k) + dt * &
973 (u_bc_accel(i,j,k) + cs%u_accel_bt(i,j,k)))
974 enddo ; enddo
975 do j=jsq,jeq ; do i=is,ie
976 v_inst(i,j,k) = g%mask2dCv(i,j) * (v_inst(i,j,k) + dt * &
977 (v_bc_accel(i,j,k) + cs%v_accel_bt(i,j,k)))
978 enddo ; enddo
979 enddo
980 call cpu_clock_end(id_clock_mom_update)
981
982 if (cs%debug) then
983 call mom_accel_chksum("Corrector accel", cs%CAu, cs%CAv, cs%PFu, cs%PFv, &
984 cs%diffu, cs%diffv, g, gv, us, cs%pbce, cs%u_accel_bt, cs%v_accel_bt, &
985 symmetric=sym)
986 call uvchksum("Corrector 1 [uv]", u_inst, v_inst, g%HI, haloshift=0, symmetric=sym, unscale=us%L_T_to_m_s)
987 call hchksum(h, "Corrector 1 h", g%HI, haloshift=1, unscale=gv%H_to_MKS)
988 call uvchksum("Corrector 1 [uv]h", uh, vh, g%HI, haloshift=2, &
989 symmetric=sym, unscale=gv%H_to_MKS*us%L_to_m**2*us%s_to_T)
990 ! call MOM_state_chksum("Corrector 1", u_inst, v_inst, h, uh, vh, G, GV, US, haloshift=1)
991
992 endif
993
994 ! u <- u + dt d/dz visc d/dz u
995 ! u_av <- u_av + dt d/dz visc d/dz u_av
996 call cpu_clock_begin(id_clock_vertvisc)
997
998 if (cs%fpmix) then
999 uold(:,:,:) = 0.0
1000 vold(:,:,:) = 0.0
1001 do k = 1, nz
1002 do j = js , je
1003 do i = isq, ieq
1004 uold(i,j,k) = u_inst(i,j,k)
1005 enddo
1006 enddo
1007 do j = jsq, jeq
1008 do i = is, ie
1009 vold(i,j,k) = v_inst(i,j,k)
1010 enddo
1011 enddo
1012 enddo
1013 endif
1014
1015 call thickness_to_dz(h, tv, dz, g, gv, us, halo_size=1)
1016 call vertvisc_coef(u_inst, v_inst, h, dz, forces, visc, tv, dt, g, gv, us, cs%vertvisc_CSp, cs%OBC, varmix)
1017
1018 if (cs%fpmix) then
1019 lfppost = .true.
1020 call vertfpmix(u_inst, v_inst, uold, vold, hbl, h, forces, dt, lfppost, cs%Cemp_NL, &
1021 g, gv, us, cs%vertvisc_CSp, cs%OBC, waves=waves)
1022 call vertvisc(u_inst, v_inst, h, forces, visc, dt, cs%OBC, cs%ADp, cs%CDp, g, gv, us, &
1023 cs%vertvisc_CSp, cs%taux_bot, cs%tauy_bot, fpmix=cs%fpmix, waves=waves)
1024
1025 else
1026 call vertvisc(u_inst, v_inst, h, forces, visc, dt, cs%OBC, cs%ADp, cs%CDp, g, gv, us, &
1027 cs%vertvisc_CSp, cs%taux_bot, cs%tauy_bot, waves=waves)
1028 endif
1029
1030 if (g%nonblocking_updates) then
1031 call cpu_clock_end(id_clock_vertvisc)
1032 call start_group_pass(cs%pass_uv, g%Domain, clock=id_clock_pass)
1033 call cpu_clock_begin(id_clock_vertvisc)
1034 endif
1035 call vertvisc_remnant(visc, cs%visc_rem_u, cs%visc_rem_v, dt, g, gv, us, cs%vertvisc_CSp)
1036 call cpu_clock_end(id_clock_vertvisc)
1037 if (showcalltree) call calltree_waypoint("done with vertvisc (step_MOM_dyn_split_RK2)")
1038
1039! Later, h_av = (h_in + h_out)/2, but for now use h_av to store h_in.
1040 !$OMP parallel do default(shared)
1041 do k=1,nz ; do j=js-cor_stencil,je+cor_stencil ; do i=is-cor_stencil,ie+cor_stencil
1042 h_av(i,j,k) = h(i,j,k)
1043 enddo ; enddo ; enddo
1044
1045 call do_group_pass(cs%pass_visc_rem, g%Domain, clock=id_clock_pass)
1046 if (g%nonblocking_updates) then
1047 call complete_group_pass(cs%pass_uv, g%Domain, clock=id_clock_pass)
1048 else
1049 call do_group_pass(cs%pass_uv, g%Domain, clock=id_clock_pass)
1050 endif
1051
1052 ! uh = u_av * h
1053 ! h = h + dt * div . uh
1054 ! u_av and v_av adjusted so their mass transports match uhbt and vhbt.
1055 call cpu_clock_begin(id_clock_continuity)
1056 call continuity(u_inst, v_inst, h, h, uh, vh, dt, g, gv, us, cs%continuity_CSp, cs%OBC, pbv, &
1057 uhbt=cs%uhbt, vhbt=cs%vhbt, visc_rem_u=cs%visc_rem_u, visc_rem_v=cs%visc_rem_v, &
1058 u_cor=u_av, v_cor=v_av)
1059 call cpu_clock_end(id_clock_continuity)
1060 call do_group_pass(cs%pass_h, g%Domain, clock=id_clock_pass)
1061 ! Whenever thickness changes let the diag manager know, target grids
1062 ! for vertical remapping may need to be regenerated.
1063 call diag_update_remap_grids(cs%diag)
1064 if (showcalltree) call calltree_waypoint("done with continuity (step_MOM_dyn_split_RK2)")
1065
1066 if (g%nonblocking_updates) then
1067 call start_group_pass(cs%pass_av_uvh, g%Domain, clock=id_clock_pass)
1068 else
1069 call do_group_pass(cs%pass_av_uvh, g%domain, clock=id_clock_pass)
1070 endif
1071
1072 if (associated(cs%OBC)) then
1073 !### I suspect that there is a bug here when u_inst is compared with a previous value of u_av
1074 ! to estimate the dominant outward group velocity, but a fix is not available yet.
1075 call radiation_open_bdry_conds(cs%OBC, u_inst, u_old_rad_obc, v_inst, v_old_rad_obc, g, gv, us, dt)
1076 endif
1077
1078! h_av = (h_in + h_out)/2 . Going in to this line, h_av = h_in.
1079 !$OMP parallel do default(shared)
1080 do k=1,nz ; do j=js-cor_stencil,je+cor_stencil ; do i=is-cor_stencil,ie+cor_stencil
1081 h_av(i,j,k) = 0.5*(h_av(i,j,k) + h(i,j,k))
1082 enddo ; enddo ; enddo
1083
1084 if (g%nonblocking_updates) &
1085 call complete_group_pass(cs%pass_av_uvh, g%Domain, clock=id_clock_pass)
1086
1087 !$OMP parallel do default(shared)
1088 do k=1,nz
1089 do j=js-2,je+2 ; do i=isq-2,ieq+2
1090 uhtr(i,j,k) = uhtr(i,j,k) + uh(i,j,k)*dt
1091 enddo ; enddo
1092 do j=jsq-2,jeq+2 ; do i=is-2,ie+2
1093 vhtr(i,j,k) = vhtr(i,j,k) + vh(i,j,k)*dt
1094 enddo ; enddo
1095 enddo
1096
1097 if (associated(cs%OBC)) then
1098 call update_segment_thickness_reservoirs(g, gv, uhtr, vhtr, h, cs%OBC)
1099 endif
1100
1101 if (cs%store_CAu) then
1102 ! Calculate a predictor-step estimate of the Coriolis and momentum advection terms
1103 ! for use in the next time step, possibly after it has been vertically remapped.
1104 call cpu_clock_begin(id_clock_cor)
1105 call disable_averaging(cs%diag) ! These calculations should not be used for diagnostics.
1106 ! CAu = -(f+zeta_av)/h_av vh + d/dx KE_av
1107 call coradcalc(u_av, v_av, h_av, uh, vh, cs%CAu_pred, cs%CAv_pred, cs%OBC, cs%AD_pred, &
1108 g, gv, us, cs%CoriolisAdv, pbv, waves=waves)
1109 cs%CAu_pred_stored = .true.
1110 call enable_averages(dt, time_local, cs%diag) ! Reenable the averaging
1111 call cpu_clock_end(id_clock_cor)
1112 if (showcalltree) call calltree_waypoint("done with CorAdCalc (step_MOM_dyn_split_RK2)")
1113 else
1114 cs%CAu_pred_stored = .false.
1115 endif
1116
1117 ! The time-averaged free surface height has already been set by the last call to btstep.
1118
1119 ! Deallocate this memory to avoid a memory leak. ### We should revisit how this array is declared. -RWH
1120 if (dyn_p_surf .and. associated(eta_pf_start)) deallocate(eta_pf_start)
1121
1122 ! Here various terms used in to update the momentum equations are
1123 ! offered for time averaging.
1124 if (cs%id_PFu > 0) call post_data(cs%id_PFu, cs%PFu, cs%diag)
1125 if (cs%id_PFv > 0) call post_data(cs%id_PFv, cs%PFv, cs%diag)
1126 if (cs%id_CAu > 0) call post_data(cs%id_CAu, cs%CAu, cs%diag)
1127 if (cs%id_CAv > 0) call post_data(cs%id_CAv, cs%CAv, cs%diag)
1128
1129 ! Here the thickness fluxes are offered for time averaging.
1130 if (cs%id_uh > 0) call post_data(cs%id_uh, uh, cs%diag)
1131 if (cs%id_vh > 0) call post_data(cs%id_vh, vh, cs%diag)
1132 if (cs%id_uav > 0) call post_data(cs%id_uav, u_av, cs%diag)
1133 if (cs%id_vav > 0) call post_data(cs%id_vav, v_av, cs%diag)
1134 if (cs%id_u_BT_accel > 0) call post_data(cs%id_u_BT_accel, cs%u_accel_bt, cs%diag)
1135 if (cs%id_v_BT_accel > 0) call post_data(cs%id_v_BT_accel, cs%v_accel_bt, cs%diag)
1136
1137 ! Calculate effective areas and post data
1138 if (cs%id_ueffA > 0) then
1139 ueffa(:,:,:) = 0
1140 do k=1,nz ; do j=js,je ; do i=isq,ieq
1141 if (abs(up(i,j,k)) > 0.) ueffa(i,j,k) = uh(i,j,k) / up(i,j,k)
1142 enddo ; enddo ; enddo
1143 call post_data(cs%id_ueffA, ueffa, cs%diag)
1144 endif
1145
1146 if (cs%id_veffA > 0) then
1147 veffa(:,:,:) = 0
1148 do k=1,nz ; do j=jsq,jeq ; do i=is,ie
1149 if (abs(vp(i,j,k)) > 0.) veffa(i,j,k) = vh(i,j,k) / vp(i,j,k)
1150 enddo ; enddo ; enddo
1151 call post_data(cs%id_veffA, veffa, cs%diag)
1152 endif
1153
1154 ! Diagnostics of the fractional thicknesses times momentum budget terms
1155 ! 3D diagnostics hf_PFu etc. are commented because there is no clarity on proper remapping grid option.
1156 ! The code is retained for debugging purposes in the future.
1157 !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)
1158 !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)
1159 !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)
1160 !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)
1161 !if (CS%id_hf_u_BT_accel > 0) &
1162 ! call post_product_u(CS%id_hf_u_BT_accel, CS%u_accel_bt, CS%ADp%diag_hfrac_u, G, nz, CS%diag)
1163 !if (CS%id_hf_v_BT_accel > 0) &
1164 ! call post_product_v(CS%id_hf_v_BT_accel, CS%v_accel_bt, CS%ADp%diag_hfrac_v, G, nz, CS%diag)
1165
1166 ! Diagnostics for the vertical sum of layer thickness x prssure force accelerations
1167 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)
1168 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)
1169
1170 ! Diagnostics for thickness-weighted vertically averaged prssure force accelerations
1171 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)
1172 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)
1173
1174 ! Diagnostics for thickness x prssure force accelerations
1175 if (cs%id_h_PFu > 0) call post_product_u(cs%id_h_PFu, cs%PFu, cs%ADp%diag_hu, g, nz, cs%diag)
1176 if (cs%id_h_PFv > 0) call post_product_v(cs%id_h_PFv, cs%PFv, cs%ADp%diag_hv, g, nz, cs%diag)
1177
1178 ! Diagnostics of Coriolis acceleratations
1179 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)
1180 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)
1181 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)
1182 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)
1183 if (cs%id_h_CAu > 0) call post_product_u(cs%id_h_CAu, cs%CAu, cs%ADp%diag_hu, g, nz, cs%diag)
1184 if (cs%id_h_CAv > 0) call post_product_v(cs%id_h_CAv, cs%CAv, cs%ADp%diag_hv, g, nz, cs%diag)
1185
1186 ! Diagnostics of barotropic solver acceleratations
1187 if (cs%id_intz_u_BT_accel_2d > 0) &
1188 call post_product_sum_u(cs%id_intz_u_BT_accel_2d, cs%u_accel_bt, cs%ADp%diag_hu, g, nz, cs%diag)
1189 if (cs%id_intz_v_BT_accel_2d > 0) &
1190 call post_product_sum_v(cs%id_intz_v_BT_accel_2d, cs%v_accel_bt, cs%ADp%diag_hv, g, nz, cs%diag)
1191 if (cs%id_hf_u_BT_accel_2d > 0) &
1192 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)
1193 if (cs%id_hf_v_BT_accel_2d > 0) &
1194 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)
1195 if (cs%id_h_u_BT_accel > 0) &
1196 call post_product_u(cs%id_h_u_BT_accel, cs%u_accel_bt, cs%ADp%diag_hu, g, nz, cs%diag)
1197 if (cs%id_h_v_BT_accel > 0) &
1198 call post_product_v(cs%id_h_v_BT_accel, cs%v_accel_bt, cs%ADp%diag_hv, g, nz, cs%diag)
1199
1200 ! Diagnostics for momentum budget terms multiplied by visc_rem_[uv],
1201 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)
1202 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)
1203 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)
1204 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)
1205 if (cs%id_u_BT_accel_visc_rem > 0) &
1206 call post_product_u(cs%id_u_BT_accel_visc_rem, cs%u_accel_bt, cs%ADp%visc_rem_u, g, nz, cs%diag)
1207 if (cs%id_v_BT_accel_visc_rem > 0) &
1208 call post_product_v(cs%id_v_BT_accel_visc_rem, cs%v_accel_bt, cs%ADp%visc_rem_v, g, nz, cs%diag)
1209
1210 ! Diagnostics related to changes in eta
1211 if (cs%id_deta_dt > 0) call post_data(cs%id_deta_dt, deta_dt, cs%diag)
1212
1213 if (cs%debug) then
1214 call mom_state_chksum("Corrector ", u_inst, v_inst, h, uh, vh, g, gv, us, symmetric=sym)
1215 call uvchksum("Corrector avg [uv]", u_av, v_av, g%HI, haloshift=1, symmetric=sym, unscale=us%L_T_to_m_s)
1216 call hchksum(h_av, "Corrector avg h", g%HI, haloshift=1, unscale=gv%H_to_MKS)
1217 ! call MOM_state_chksum("Corrector avg ", u_av, v_av, h_av, uh, vh, G, GV, US)
1218 endif
1219
1220 if (showcalltree) call calltree_leave("step_MOM_dyn_split_RK2()")
1221
1222end subroutine step_mom_dyn_split_rk2
1223
1224!> This subroutine sets up any auxiliary restart variables that are specific
1225!! to the split-explicit time stepping scheme. All variables registered here should
1226!! have the ability to be recreated if they are not present in a restart file.
1227subroutine register_restarts_dyn_split_rk2(HI, GV, US, param_file, CS, restart_CS, uh, vh)
1228 type(hor_index_type), intent(in) :: hi !< Horizontal index structure
1229 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
1230 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1231 type(param_file_type), intent(in) :: param_file !< parameter file
1232 type(mom_dyn_split_rk2_cs), pointer :: cs !< module control structure
1233 type(mom_restart_cs), intent(inout) :: restart_cs !< MOM restart control structure
1234 real, dimension(SZIB_(HI),SZJ_(HI),SZK_(GV)), &
1235 target, intent(inout) :: uh !< zonal volume or mass transport [H L2 T-1 ~> m3 s-1 or kg s-1]
1236 real, dimension(SZI_(HI),SZJB_(HI),SZK_(GV)), &
1237 target, intent(inout) :: vh !< merid volume or mass transport [H L2 T-1 ~> m3 s-1 or kg s-1]
1238
1239 character(len=40) :: mdl = "MOM_dynamics_split_RK2" ! This module's name.
1240 type(vardesc) :: vd(2)
1241 character(len=48) :: thickness_units, flux_units
1242
1243 integer :: isd, ied, jsd, jed, nz, isdb, iedb, jsdb, jedb
1244
1245 isd = hi%isd ; ied = hi%ied ; jsd = hi%jsd ; jed = hi%jed ; nz = gv%ke
1246 isdb = hi%IsdB ; iedb = hi%IedB ; jsdb = hi%JsdB ; jedb = hi%JedB
1247
1248 ! This is where a control structure specific to this module would be allocated.
1249 if (associated(cs)) then
1250 call mom_error(warning, "register_restarts_dyn_split_RK2 called with an associated "// &
1251 "control structure.")
1252 return
1253 endif
1254 allocate(cs)
1255
1256 alloc_(cs%diffu(isdb:iedb,jsd:jed,nz)) ; cs%diffu(:,:,:) = 0.0
1257 alloc_(cs%diffv(isd:ied,jsdb:jedb,nz)) ; cs%diffv(:,:,:) = 0.0
1258 alloc_(cs%CAu(isdb:iedb,jsd:jed,nz)) ; cs%CAu(:,:,:) = 0.0
1259 alloc_(cs%CAv(isd:ied,jsdb:jedb,nz)) ; cs%CAv(:,:,:) = 0.0
1260 alloc_(cs%CAu_pred(isdb:iedb,jsd:jed,nz)) ; cs%CAu_pred(:,:,:) = 0.0
1261 alloc_(cs%CAv_pred(isd:ied,jsdb:jedb,nz)) ; cs%CAv_pred(:,:,:) = 0.0
1262 alloc_(cs%PFu(isdb:iedb,jsd:jed,nz)) ; cs%PFu(:,:,:) = 0.0
1263 alloc_(cs%PFv(isd:ied,jsdb:jedb,nz)) ; cs%PFv(:,:,:) = 0.0
1264
1265 alloc_(cs%eta(isd:ied,jsd:jed)) ; cs%eta(:,:) = 0.0
1266 alloc_(cs%u_av(isdb:iedb,jsd:jed,nz)) ; cs%u_av(:,:,:) = 0.0
1267 alloc_(cs%v_av(isd:ied,jsdb:jedb,nz)) ; cs%v_av(:,:,:) = 0.0
1268 alloc_(cs%h_av(isd:ied,jsd:jed,nz)) ; cs%h_av(:,:,:) = gv%Angstrom_H
1269
1270 allocate(cs%taux_bot(isdb:iedb,jsd:jed), source = 0.0)
1271 allocate(cs%tauy_bot(isd:ied,jsdb:jedb), source = 0.0)
1272
1273 thickness_units = get_thickness_units(gv)
1274 flux_units = get_flux_units(gv)
1275
1276 call get_param(param_file, mdl, "STORE_CORIOLIS_ACCEL", cs%store_CAu, &
1277 "If true, calculate the Coriolis accelerations at the end of each "//&
1278 "timestep for use in the predictor step of the next split RK2 timestep.", &
1279 default=.true., do_not_log=.true.)
1280
1281 if (gv%Boussinesq) then
1282 call register_restart_field(cs%eta, "sfc", .false., restart_cs, &
1283 longname="Free surface Height", units=thickness_units, conversion=gv%H_to_mks)
1284 else
1285 call register_restart_field(cs%eta, "p_bot", .false., restart_cs, &
1286 longname="Bottom Pressure", units=thickness_units, conversion=gv%H_to_mks)
1287 endif
1288
1289 ! These are needed, either to calculate CAu and CAv or to calculate the velocity anomalies in
1290 ! the barotropic solver's Coriolis terms.
1291 vd(1) = var_desc("u2", "m s-1", "Auxiliary Zonal velocity", 'u', 'L')
1292 vd(2) = var_desc("v2", "m s-1", "Auxiliary Meridional velocity", 'v', 'L')
1293 call register_restart_pair(cs%u_av, cs%v_av, vd(1), vd(2), .false., restart_cs, &
1294 conversion=us%L_T_to_m_s)
1295
1296 if (cs%store_CAu) then
1297 vd(1) = var_desc("CAu", "m s-2", "Zonal Coriolis and advactive acceleration", 'u', 'L')
1298 vd(2) = var_desc("CAv", "m s-2", "Meridional Coriolis and advactive acceleration", 'v', 'L')
1299 call register_restart_pair(cs%CAu_pred, cs%CAv_pred, vd(1), vd(2), .false., restart_cs, &
1300 conversion=us%L_T2_to_m_s2)
1301 else
1302 call register_restart_field(cs%h_av, "h2", .false., restart_cs, &
1303 longname="Auxiliary Layer Thickness", units=thickness_units, conversion=gv%H_to_mks)
1304
1305 vd(1) = var_desc("uh", flux_units, "Zonal thickness flux", 'u', 'L')
1306 vd(2) = var_desc("vh", flux_units, "Meridional thickness flux", 'v', 'L')
1307 call register_restart_pair(uh, vh, vd(1), vd(2), .false., restart_cs, &
1308 conversion=gv%H_to_MKS*us%L_to_m**2*us%s_to_T)
1309 endif
1310
1311 vd(1) = var_desc("diffu", "m s-2", "Zonal horizontal viscous acceleration", 'u', 'L')
1312 vd(2) = var_desc("diffv", "m s-2", "Meridional horizontal viscous acceleration", 'v', 'L')
1313 call register_restart_pair(cs%diffu, cs%diffv, vd(1), vd(2), .false., restart_cs, &
1314 conversion=us%L_T2_to_m_s2)
1315
1316 call register_barotropic_restarts(hi, gv, us, param_file, cs%barotropic_CSp, restart_cs)
1317
1318 call get_param(param_file, mdl, "SPLIT_BOTTOM_STRESS", cs%split_bottom_stress, &
1319 "If true, provide the bottom stress calculated by the "//&
1320 "vertical viscosity to the barotropic solver.", default=.false.,&
1321 do_not_log=.true.)
1322
1323 if (cs%split_bottom_stress) then
1324 vd(1) = var_desc("taux_bot", "kg m-1 s-2", "Zonal bottom stress", 'u', '1')
1325 vd(2) = var_desc("tauy_bot", "kg m-1 s-2", "Meridional bottom stress", 'v', '1')
1326 call register_restart_pair(cs%taux_bot, cs%tauy_bot, vd(1), vd(2), .false., restart_cs, &
1327 conversion=us%RLZ_T2_to_Pa)
1328 endif
1329
1331
1332!> This subroutine does remapping for the auxiliary restart variables that are used
1333!! with the split RK2 time stepping scheme.
1334subroutine remap_dyn_split_rk2_aux_vars(G, GV, CS, h_old_u, h_old_v, h_new_u, h_new_v, ALE_CSp)
1335 type(ocean_grid_type), intent(inout) :: g !< ocean grid structure
1336 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
1337 type(mom_dyn_split_rk2_cs), pointer :: cs !< module control structure
1338 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
1339 intent(in) :: h_old_u !< Source grid thickness at zonal
1340 !! velocity points [H ~> m or kg m-2]
1341 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
1342 intent(in) :: h_old_v !< Source grid thickness at meridional
1343 !! velocity points [H ~> m or kg m-2]
1344 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
1345 intent(in) :: h_new_u !< Destination grid thickness at zonal
1346 !! velocity points [H ~> m or kg m-2]
1347 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
1348 intent(in) :: h_new_v !< Destination grid thickness at meridional
1349 !! velocity points [H ~> m or kg m-2]
1350 type(ale_cs), pointer :: ale_csp !< ALE control structure to use when remapping
1351
1352 if (.not.cs%remap_aux) return
1353
1354 if (cs%store_CAu) then
1355 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)
1356 call pass_vector(cs%u_av, cs%v_av, g%Domain, complete=.false.)
1357 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)
1358 call pass_vector(cs%CAu_pred, cs%CAv_pred, g%Domain, complete=.false.)
1359 endif
1360
1361 call ale_remap_velocities(ale_csp, g, gv, h_old_u, h_old_v, h_new_u, h_new_v, cs%diffu, cs%diffv)
1362 call pass_vector(cs%diffu, cs%diffv, g%Domain, complete=.true.)
1363
1364end subroutine remap_dyn_split_rk2_aux_vars
1365
1366!> Initializes aspects of the dyn_split_RK2 that depend on diabatic processes.
1367!! Needed when BLDs are used in the dynamics.
1368subroutine init_dyn_split_rk2_diabatic(diabatic_CSp, CS)
1369 type(diabatic_cs), intent(in) :: diabatic_csp !< diabatic structure
1370 type(mom_dyn_split_rk2_cs), pointer :: cs !< module control structure
1371
1372 call extract_diabatic_member(diabatic_csp, kpp_csp=cs%KPP_CSp)
1373 call extract_diabatic_member(diabatic_csp, energetic_pbl_csp=cs%energetic_PBL_CSp)
1374
1375end subroutine init_dyn_split_rk2_diabatic
1376
1377!> This subroutine initializes all of the variables that are used by this
1378!! dynamic core, including diagnostics and the cpu clocks.
1379subroutine initialize_dyn_split_rk2(u, v, h, tv, uh, vh, eta, Time, G, GV, US, param_file, &
1380 diag, CS, HA_CSp, restart_CS, dt, Accel_diag, Cont_diag, MIS, &
1381 VarMix, MEKE, thickness_diffuse_CSp, &
1382 OBC, update_OBC_CSp, ALE_CSp, set_visc, &
1383 visc, dirs, ntrunc, pbv, calc_dtbt, cont_stencil, dyn_h_stencil)
1384 type(ocean_grid_type), intent(inout) :: g !< ocean grid structure
1385 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
1386 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1387 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
1388 intent(inout) :: u !< zonal velocity [L T-1 ~> m s-1]
1389 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
1390 intent(inout) :: v !< merid velocity [L T-1 ~> m s-1]
1391 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1392 intent(inout) :: h !< layer thickness [H ~> m or kg m-2]
1393 type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic type
1394 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
1395 target, intent(inout) :: uh !< zonal volume/mass transport [H L2 T-1 ~> m3 s-1 or kg s-1]
1396 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
1397 target, intent(inout) :: vh !< merid volume/mass transport [H L2 T-1 ~> m3 s-1 or kg s-1]
1398 real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: eta !< free surface height or column mass [H ~> m or kg m-2]
1399 type(time_type), target, intent(in) :: time !< current model time
1400 type(param_file_type), intent(in) :: param_file !< parameter file for parsing
1401 type(diag_ctrl), target, intent(inout) :: diag !< to control diagnostics
1402 type(mom_dyn_split_rk2_cs), pointer :: cs !< module control structure
1403 type(harmonic_analysis_cs), pointer :: ha_csp !< A pointer to the control structure of the
1404 !! harmonic analysis module
1405 type(mom_restart_cs), intent(inout) :: restart_cs !< MOM restart control structure
1406 real, intent(in) :: dt !< time step [T ~> s]
1407 type(accel_diag_ptrs), target, intent(inout) :: accel_diag !< points to momentum equation terms for
1408 !! budget analysis
1409 type(cont_diag_ptrs), target, intent(inout) :: cont_diag !< points to terms in continuity equation
1410 type(ocean_internal_state), intent(inout) :: mis !< "MOM6 internal state" used to pass
1411 !! diagnostic pointers
1412 type(varmix_cs), intent(inout) :: varmix !< points to spatially variable viscosities
1413 type(meke_type), intent(inout) :: meke !< MEKE fields
1414 type(thickness_diffuse_cs), intent(inout) :: thickness_diffuse_csp !< Pointer to the control structure
1415 !! used for the isopycnal height diffusive transport.
1416 type(ocean_obc_type), pointer :: obc !< points to OBC related fields
1417 type(update_obc_cs), pointer :: update_obc_csp !< points to OBC update related fields
1418 type(ale_cs), pointer :: ale_csp !< points to ALE control structure
1419 type(set_visc_cs), target, intent(in) :: set_visc !< set_visc control structure
1420 type(vertvisc_type), intent(inout) :: visc !< vertical viscosities, bottom drag, and related
1421 type(directories), intent(in) :: dirs !< contains directory paths
1422 integer, target, intent(inout) :: ntrunc !< A target for the variable that records
1423 !! the number of times the velocity is
1424 !! truncated (this should be 0).
1425 logical, intent(out) :: calc_dtbt !< If true, recalculate the barotropic time step
1426 type(porous_barrier_type), intent(in) :: pbv !< porous barrier fractional cell metrics
1427 integer, intent(out) :: cont_stencil !< The stencil for thickness
1428 !! from the continuity solver.
1429 integer, intent(out) :: dyn_h_stencil !< The stencil for thickness for the
1430 !! the dynamics based on the continuity
1431 !! solver and Coriolis scheme.
1432
1433 ! local variables
1434 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_tmp ! A temporary copy of the layer thicknesses [H ~> m or kg m-2]
1435 character(len=40) :: mdl = "MOM_dynamics_split_RK2" ! This module's name.
1436 ! This include declares and sets the variable "version".
1437# include "version_variable.h"
1438 character(len=48) :: thickness_units, flux_units, eta_rest_name
1439 type(group_pass_type) :: pass_av_h_uvh
1440 logical :: debug_truncations
1441 logical :: read_uv, read_h2
1442 logical :: enable_bugs ! If true, the defaults for recently added bug-fix flags are set to
1443 ! recreate the bugs, or if false bugs are only used if actively selected.
1444 logical :: visc_rem_bug ! Stores the value of runtime paramter VISC_REM_BUG.
1445 integer :: cor_stencil
1446
1447 integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz
1448 integer :: isdb, iedb, jsdb, jedb
1449 integer :: nc ! Number of tidal constituents to be harmonically analyzed
1450 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
1451 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
1452 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
1453
1454 if (.not.associated(cs)) call mom_error(fatal, &
1455 "initialize_dyn_split_RK2 called with an unassociated control structure.")
1456 if (cs%module_is_initialized) then
1457 call mom_error(warning, "initialize_dyn_split_RK2 called with a control "// &
1458 "structure that has already been initialized.")
1459 return
1460 endif
1461 cs%module_is_initialized = .true.
1462
1463 cs%diag => diag
1464
1465 call log_version(param_file, mdl, version, "")
1466 call get_param(param_file, mdl, "TIDES", cs%use_tides, &
1467 "If true, apply tidal momentum forcing.", default=.false.)
1468 call get_param(param_file, mdl, "CALCULATE_SAL", cs%calculate_SAL, &
1469 "If true, calculate self-attraction and loading.", default=cs%use_tides)
1470 call get_param(param_file, mdl, "USE_HA", cs%use_HA, &
1471 "If true, perform inline harmonic analysis.", default=.false.)
1472 call get_param(param_file, mdl, "HA_N_CONST", nc, &
1473 "Number of tidal constituents to be harmonically analyzed.", &
1474 default=0, do_not_log=.not.cs%use_HA)
1475 if (nc<=0) cs%use_HA = .false.
1476 call get_param(param_file, mdl, "BE", cs%be, &
1477 "If SPLIT is true, BE determines the relative weighting "//&
1478 "of a 2nd-order Runga-Kutta baroclinic time stepping "//&
1479 "scheme (0.5) and a backward Euler scheme (1) that is "//&
1480 "used for the Coriolis and inertial terms. BE may be "//&
1481 "from 0.5 to 1, but instability may occur near 0.5. "//&
1482 "BE is also applicable if SPLIT is false and USE_RK2 "//&
1483 "is true.", units="nondim", default=0.6)
1484 call get_param(param_file, mdl, "BEGW", cs%begw, &
1485 "If SPLIT is true, BEGW is a number from 0 to 1 that "//&
1486 "controls the extent to which the treatment of gravity "//&
1487 "waves is forward-backward (0) or simulated backward "//&
1488 "Euler (1). 0 is almost always used. "//&
1489 "If SPLIT is false and USE_RK2 is true, BEGW can be "//&
1490 "between 0 and 0.5 to damp gravity waves.", &
1491 units="nondim", default=0.0)
1492 call get_param(param_file, mdl, "SET_DTBT_USE_BT_CONT", cs%dtbt_use_bt_cont, &
1493 "If true, use BT_CONT to calculate DTBT if possible.", default=.false.)
1494 call get_param(param_file, mdl, "SPLIT_BOTTOM_STRESS", cs%split_bottom_stress, &
1495 "If true, provide the bottom stress calculated by the "//&
1496 "vertical viscosity to the barotropic solver.", default=.false.)
1497 call get_param(param_file, mdl, "BT_USE_LAYER_FLUXES", cs%BT_use_layer_fluxes, &
1498 "If true, use the summed layered fluxes plus an "//&
1499 "adjustment due to the change in the barotropic velocity "//&
1500 "in the barotropic continuity equation.", default=.true.)
1501 call get_param(param_file, mdl, "BT_ADJ_CORR_MASS_SRC", cs%BT_adj_corr_mass_src, &
1502 "If true, recalculates the barotropic mass source after "//&
1503 "predictor step. This should make little difference in the "//&
1504 "deep ocean but appears to help for vanished layers. If false, "//&
1505 "uses the same mass source as from the predictor step.", default=.true.)
1506 call get_param(param_file, mdl, "STORE_CORIOLIS_ACCEL", cs%store_CAu, &
1507 "If true, calculate the Coriolis accelerations at the end of each "//&
1508 "timestep for use in the predictor step of the next split RK2 timestep.", &
1509 default=.true.)
1510 call get_param(param_file, mdl, "FPMIX", cs%fpmix, &
1511 "If true, add non-local momentum flux increments and diffuse down the Eulerian gradient.", &
1512 default=.false.)
1513 if (cs%fpmix) then
1514 call get_param(param_file, "MOM", "CEMP_NL", cs%Cemp_NL, &
1515 "Empirical coefficient of non-local momentum mixing.", &
1516 units="nondim", default=3.6)
1517 endif
1518 call get_param(param_file, mdl, "REMAP_AUXILIARY_VARS", cs%remap_aux, &
1519 "If true, apply ALE remapping to all of the auxiliary 3-dimensional "//&
1520 "variables that are needed to reproduce across restarts, similarly to "//&
1521 "what is already being done with the primary state variables. "//&
1522 "The default should be changed to true.", default=.false., do_not_log=.true.)
1523 if (cs%remap_aux .and. .not.cs%store_CAu) call mom_error(fatal, &
1524 "REMAP_AUXILIARY_VARS requires that STORE_CORIOLIS_ACCEL = True.")
1525 call get_param(param_file, mdl, "DEBUG", cs%debug, &
1526 "If true, write out verbose debugging data.", &
1527 default=.false., debuggingparam=.true.)
1528 call get_param(param_file, mdl, "OBC_DEBUGGING_TESTS", cs%debug_OBC, &
1529 "If true, do additional calls resetting certain values to help verify the "//&
1530 "correctness of the open boundary condition code.", &
1531 default=.false., old_name="DEBUG_OBC", debuggingparam=.true., do_not_log=.true.)
1532 call get_param(param_file, mdl, "DEBUG_TRUNCATIONS", debug_truncations, &
1533 default=.false.)
1534 call get_param(param_file, mdl, "ENABLE_BUGS_BY_DEFAULT", enable_bugs, &
1535 default=.true., do_not_log=.true.) ! This is logged from MOM.F90.
1536 call get_param(param_file, mdl, "VISC_REM_BUG", visc_rem_bug, &
1537 "If true, visc_rem_[uv] in split mode is incorrectly calculated or accounted "//&
1538 "for in two places. This parameter controls the defaults of two individual "//&
1539 "flags, VISC_REM_TIMESTEP_BUG in MOM_dynamics_split_RK2(b) and "//&
1540 "VISC_REM_BT_WEIGHT_BUG in MOM_barotropic.", default=.false.)
1541 call get_param(param_file, mdl, "VISC_REM_TIMESTEP_BUG", cs%visc_rem_dt_bug, &
1542 "If true, recover a bug that uses dt_pred rather than dt in "//&
1543 "vertvisc_remnant() at the end of predictor stage for the following "//&
1544 "continuity() and btstep() calls in the corrector step. Default of this flag "//&
1545 "is set by VISC_REM_BUG", default=visc_rem_bug)
1546
1547 alloc_(cs%uhbt(isdb:iedb,jsd:jed)) ; cs%uhbt(:,:) = 0.0
1548 alloc_(cs%vhbt(isd:ied,jsdb:jedb)) ; cs%vhbt(:,:) = 0.0
1549 alloc_(cs%visc_rem_u(isdb:iedb,jsd:jed,nz)) ; cs%visc_rem_u(:,:,:) = 0.0
1550 alloc_(cs%visc_rem_v(isd:ied,jsdb:jedb,nz)) ; cs%visc_rem_v(:,:,:) = 0.0
1551 alloc_(cs%eta_PF(isd:ied,jsd:jed)) ; cs%eta_PF(:,:) = 0.0
1552 alloc_(cs%pbce(isd:ied,jsd:jed,nz)) ; cs%pbce(:,:,:) = 0.0
1553
1554 alloc_(cs%u_accel_bt(isdb:iedb,jsd:jed,nz)) ; cs%u_accel_bt(:,:,:) = 0.0
1555 alloc_(cs%v_accel_bt(isd:ied,jsdb:jedb,nz)) ; cs%v_accel_bt(:,:,:) = 0.0
1556 alloc_(cs%PFu_Stokes(isdb:iedb,jsd:jed,nz)) ; cs%PFu_Stokes(:,:,:) = 0.0
1557 alloc_(cs%PFv_Stokes(isd:ied,jsdb:jedb,nz)) ; cs%PFv_Stokes(:,:,:) = 0.0
1558
1559 mis%diffu => cs%diffu
1560 mis%diffv => cs%diffv
1561 mis%PFu => cs%PFu
1562 mis%PFv => cs%PFv
1563 mis%CAu => cs%CAu
1564 mis%CAv => cs%CAv
1565 mis%pbce => cs%pbce
1566 mis%u_accel_bt => cs%u_accel_bt
1567 mis%v_accel_bt => cs%v_accel_bt
1568 mis%u_av => cs%u_av
1569 mis%v_av => cs%v_av
1570
1571 cs%ADp => accel_diag
1572 cs%CDp => cont_diag
1573 accel_diag%diffu => cs%diffu
1574 accel_diag%diffv => cs%diffv
1575 accel_diag%PFu => cs%PFu
1576 accel_diag%PFv => cs%PFv
1577 accel_diag%CAu => cs%CAu
1578 accel_diag%CAv => cs%CAv
1579 accel_diag%u_accel_bt => cs%u_accel_bt
1580 accel_diag%v_accel_bt => cs%v_accel_bt
1581
1582 allocate(cs%AD_pred)
1583 cs%AD_pred%diffu => cs%diffu
1584 cs%AD_pred%diffv => cs%diffv
1585 cs%AD_pred%PFu => cs%PFu
1586 cs%AD_pred%PFv => cs%PFv
1587 cs%AD_pred%CAu => cs%CAu_pred
1588 cs%AD_pred%CAv => cs%CAv_pred
1589 cs%AD_pred%u_accel_bt => cs%u_accel_bt
1590 cs%AD_pred%v_accel_bt => cs%v_accel_bt
1591
1592! Accel_diag%pbce => CS%pbce
1593! Accel_diag%u_accel_bt => CS%u_accel_bt ; Accel_diag%v_accel_bt => CS%v_accel_bt
1594! Accel_diag%u_av => CS%u_av ; Accel_diag%v_av => CS%v_av
1595
1596 id_clock_pass_init = cpu_clock_id('(Ocean init message passing)', grain=clock_routine)
1597
1598 call continuity_init(time, g, gv, us, param_file, diag, cs%continuity_CSp, cs%OBC)
1599 cont_stencil = continuity_stencil(cs%continuity_CSp)
1600 call coriolisadv_init(time, g, gv, us, param_file, diag, cs%ADp, cs%CoriolisAdv)
1601 cor_stencil = coriolisadv_stencil(cs%CoriolisAdv)
1602 dyn_h_stencil = max(cont_stencil, coriolisadv_stencil(cs%CoriolisAdv))
1603 if (cs%calculate_SAL) call sal_init(h, tv, g, gv, us, param_file, cs%SAL_CSp, restart_cs)
1604 if (cs%use_tides) call tidal_forcing_init(time, g, us, param_file, cs%tides_CSp)
1605 if (cs%use_HA) then
1606 call ha_init(time, us, param_file, nc, cs%HA_CSp)
1607 ha_csp => cs%HA_CSp
1608 else
1609 ha_csp => null()
1610 endif
1611 call pressureforce_init(time, g, gv, us, param_file, diag, cs%PressureForce_CSp, cs%ADp, &
1612 cs%SAL_CSp, cs%tides_CSp)
1613 call hor_visc_init(time, g, gv, us, param_file, diag, cs%hor_visc, adp=cs%ADp)
1614 call vertvisc_init(mis, time, g, gv, us, param_file, diag, cs%ADp, dirs, &
1615 ntrunc, cs%vertvisc_CSp, cs%fpmix)
1616 cs%set_visc_CSp => set_visc
1617 call updatecfltruncationvalue(time, cs%vertvisc_CSp, us, activate=is_new_run(restart_cs) )
1618
1619 if (associated(ale_csp)) cs%ALE_CSp => ale_csp
1620 if (associated(obc)) then
1621 cs%OBC => obc
1622 if (obc%ramp) call update_obc_ramp(time, cs%OBC, us, activate=is_new_run(restart_cs) )
1623 endif
1624 if (associated(update_obc_csp)) cs%update_OBC_CSp => update_obc_csp
1625
1626 eta_rest_name = "sfc" ; if (.not.gv%Boussinesq) eta_rest_name = "p_bot"
1627 if (.not. query_initialized(cs%eta, trim(eta_rest_name), restart_cs)) then
1628 ! Estimate eta based on the layer thicknesses - h. With the Boussinesq
1629 ! approximation, eta is the free surface height anomaly, while without it
1630 ! eta is the mass of ocean per unit area. eta always has the same
1631 ! dimensions as h, either m or kg m-3.
1632 ! CS%eta(:,:) = 0.0 already from initialization.
1633 if (gv%Boussinesq) then
1634 do j=js,je ; do i=is,ie ; cs%eta(i,j) = -gv%Z_to_H * g%bathyT(i,j) ; enddo ; enddo
1635 endif
1636 do k=1,nz ; do j=js,je ; do i=is,ie
1637 cs%eta(i,j) = cs%eta(i,j) + h(i,j,k)
1638 enddo ; enddo ; enddo
1639 call set_initialized(cs%eta, trim(eta_rest_name), restart_cs)
1640 endif
1641 ! Copy eta into an output array.
1642 do j=js,je ; do i=is,ie ; eta(i,j) = cs%eta(i,j) ; enddo ; enddo
1643
1644 call barotropic_init(u, v, h, time, g, gv, us, param_file, diag, &
1645 cs%barotropic_CSp, restart_cs, calc_dtbt, cs%BT_cont, &
1646 cs%OBC, cs%SAL_CSp, ha_csp)
1647
1648 if (.not. query_initialized(cs%diffu, "diffu", restart_cs) .or. &
1649 .not. query_initialized(cs%diffv, "diffv", restart_cs)) then
1650 call horizontal_viscosity(u, v, h, uh, vh, cs%diffu, cs%diffv, meke, varmix, g, gv, us, cs%hor_visc, &
1651 tv, dt, obc=cs%OBC, bt=cs%barotropic_CSp, td=thickness_diffuse_csp, &
1652 hu_cont=cs%BT_cont%h_u, hv_cont=cs%BT_cont%h_v)
1653 call set_initialized(cs%diffu, "diffu", restart_cs)
1654 call set_initialized(cs%diffv, "diffv", restart_cs)
1655 endif
1656
1657 if (.not. query_initialized(cs%u_av, "u2", restart_cs) .or. &
1658 .not. query_initialized(cs%v_av, "v2", restart_cs)) then
1659 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
1660 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
1661 call set_initialized(cs%u_av, "u2", restart_cs)
1662 call set_initialized(cs%v_av, "v2", restart_cs)
1663 endif
1664
1665 if (cs%store_CAu) then
1666 if (query_initialized(cs%CAu_pred, "CAu", restart_cs) .and. &
1667 query_initialized(cs%CAv_pred, "CAv", restart_cs)) then
1668 cs%CAu_pred_stored = .true.
1669 else
1670 call only_read_from_restarts(uh, vh, 'uh', 'vh', g, restart_cs, stagger=cgrid_ne, &
1671 filename=dirs%input_filename, directory=dirs%restart_input_dir, &
1672 success=read_uv, scale=us%m_to_L**2*us%T_to_s/gv%H_to_mks)
1673 call only_read_from_restarts('h2', cs%h_av, g, restart_cs, &
1674 filename=dirs%input_filename, directory=dirs%restart_input_dir, &
1675 success=read_h2, scale=1.0/gv%H_to_mks)
1676 if (read_uv .and. read_h2) then
1677 call pass_var(cs%h_av, g%Domain, clock=id_clock_pass_init)
1678 else
1679 do k=1,nz ; do j=jsd,jed ; do i=isd,ied ; h_tmp(i,j,k) = h(i,j,k) ; enddo ; enddo ; enddo
1680 call continuity(cs%u_av, cs%v_av, h, h_tmp, uh, vh, dt, g, gv, us, cs%continuity_CSp, cs%OBC, pbv)
1681 call pass_var(h_tmp, g%Domain, clock=id_clock_pass_init)
1682 do k=1,nz ; do j=jsd,jed ; do i=isd,ied
1683 cs%h_av(i,j,k) = 0.5*(h(i,j,k) + h_tmp(i,j,k))
1684 enddo ; enddo ; enddo
1685 endif
1686 call pass_vector(cs%u_av, cs%v_av, g%Domain, halo=cor_stencil, clock=id_clock_pass_init, complete=.false.)
1687 call pass_vector(uh, vh, g%Domain, halo=cor_stencil, clock=id_clock_pass_init, complete=.true.)
1688 call coradcalc(cs%u_av, cs%v_av, cs%h_av, uh, vh, cs%CAu_pred, cs%CAv_pred, cs%OBC, cs%ADp, &
1689 g, gv, us, cs%CoriolisAdv, pbv) !, Waves=Waves)
1690 cs%CAu_pred_stored = .true.
1691 endif
1692 else
1693 cs%CAu_pred_stored = .false.
1694 ! This call is just here to initialize uh and vh.
1695 if (.not. query_initialized(uh, "uh", restart_cs) .or. &
1696 .not. query_initialized(vh, "vh", restart_cs)) then
1697 do k=1,nz ; do j=jsd,jed ; do i=isd,ied ; h_tmp(i,j,k) = h(i,j,k) ; enddo ; enddo ; enddo
1698 call continuity(u, v, h, h_tmp, uh, vh, dt, g, gv, us, cs%continuity_CSp, cs%OBC, pbv)
1699 call pass_var(h_tmp, g%Domain, clock=id_clock_pass_init)
1700 do k=1,nz ; do j=jsd,jed ; do i=isd,ied
1701 cs%h_av(i,j,k) = 0.5*(h(i,j,k) + h_tmp(i,j,k))
1702 enddo ; enddo ; enddo
1703 call set_initialized(uh, "uh", restart_cs)
1704 call set_initialized(vh, "vh", restart_cs)
1705 call set_initialized(cs%h_av, "h2", restart_cs)
1706 ! Try reading the CAu and CAv fields from the restart file, in case this restart file is
1707 ! using a newer format.
1708 call only_read_from_restarts(cs%CAu_pred, cs%CAv_pred, "CAu", "CAv", g, restart_cs, &
1709 stagger=cgrid_ne, filename=dirs%input_filename, directory=dirs%restart_input_dir, &
1710 success=read_uv, scale=us%m_s_to_L_T*us%T_to_s)
1711 cs%CAu_pred_stored = read_uv
1712 else
1713 if (.not. query_initialized(cs%h_av, "h2", restart_cs)) then
1714 cs%h_av(:,:,:) = h(:,:,:)
1715 call set_initialized(cs%h_av, "h2", restart_cs)
1716 endif
1717 endif
1718 endif
1719 call cpu_clock_begin(id_clock_pass_init)
1720 call create_group_pass(pass_av_h_uvh, cs%u_av, cs%v_av, g%Domain, halo=2)
1721 if (cs%CAu_pred_stored) then
1722 call create_group_pass(pass_av_h_uvh, cs%CAu_pred, cs%CAv_pred, g%Domain, halo=2)
1723 else
1724 call create_group_pass(pass_av_h_uvh, cs%h_av, g%Domain, halo=2)
1725 call create_group_pass(pass_av_h_uvh, uh, vh, g%Domain, halo=2)
1726 endif
1727 call do_group_pass(pass_av_h_uvh, g%Domain)
1728 call cpu_clock_end(id_clock_pass_init)
1729
1730 flux_units = get_flux_units(gv)
1731 thickness_units = get_thickness_units(gv)
1732 cs%id_uh = register_diag_field('ocean_model', 'uh', diag%axesCuL, time, &
1733 'Zonal Thickness Flux', flux_units, conversion=gv%H_to_MKS*us%L_to_m**2*us%s_to_T, &
1734 y_cell_method='sum', v_extensive=.true.)
1735 cs%id_vh = register_diag_field('ocean_model', 'vh', diag%axesCvL, time, &
1736 'Meridional Thickness Flux', flux_units, conversion=gv%H_to_MKS*us%L_to_m**2*us%s_to_T, &
1737 x_cell_method='sum', v_extensive=.true.)
1738
1739 cs%id_CAu = register_diag_field('ocean_model', 'CAu', diag%axesCuL, time, &
1740 'Zonal Coriolis and Advective Acceleration', 'm s-2', conversion=us%L_T2_to_m_s2)
1741 cs%id_CAv = register_diag_field('ocean_model', 'CAv', diag%axesCvL, time, &
1742 'Meridional Coriolis and Advective Acceleration', 'm s-2', conversion=us%L_T2_to_m_s2)
1743 cs%id_PFu = register_diag_field('ocean_model', 'PFu', diag%axesCuL, time, &
1744 'Zonal Pressure Force Acceleration', 'm s-2', conversion=us%L_T2_to_m_s2)
1745 cs%id_PFv = register_diag_field('ocean_model', 'PFv', diag%axesCvL, time, &
1746 'Meridional Pressure Force Acceleration', 'm s-2', conversion=us%L_T2_to_m_s2)
1747 cs%id_ueffA = register_diag_field('ocean_model', 'ueffA', diag%axesCuL, time, &
1748 'Effective U-Face Area', 'm^2', conversion=gv%H_to_m*us%L_to_m, &
1749 y_cell_method='sum', v_extensive=.true.)
1750 cs%id_veffA = register_diag_field('ocean_model', 'veffA', diag%axesCvL, time, &
1751 'Effective V-Face Area', 'm^2', conversion=gv%H_to_m*us%L_to_m, &
1752 x_cell_method='sum', v_extensive=.true.)
1753 if (gv%Boussinesq) then
1754 cs%id_deta_dt = register_diag_field('ocean_model', 'deta_dt', diag%axesT1, time, &
1755 'Barotropic SSH tendency due to dynamics', trim(thickness_units)//' s-1', conversion=gv%H_to_MKS*us%s_to_T)
1756 else
1757 cs%id_deta_dt = register_diag_field('ocean_model', 'deta_dt', diag%axesT1, time, &
1758 'Barotropic column-mass tendency due to dynamics', trim(thickness_units)//' s-1', &
1759 conversion=gv%H_to_mks*us%s_to_T)
1760 endif
1761
1762 !CS%id_hf_PFu = register_diag_field('ocean_model', 'hf_PFu', diag%axesCuL, Time, &
1763 ! 'Fractional Thickness-weighted Zonal Pressure Force Acceleration', &
1764 ! 'm s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)
1765 !if (CS%id_hf_PFu > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_u,IsdB,IedB,jsd,jed,nz)
1766
1767 !CS%id_hf_PFv = register_diag_field('ocean_model', 'hf_PFv', diag%axesCvL, Time, &
1768 ! 'Fractional Thickness-weighted Meridional Pressure Force Acceleration', &
1769 ! 'm s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)
1770 !if (CS%id_hf_PFv > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz)
1771
1772 !CS%id_hf_CAu = register_diag_field('ocean_model', 'hf_CAu', diag%axesCuL, Time, &
1773 ! 'Fractional Thickness-weighted Zonal Coriolis and Advective Acceleration', &
1774 ! 'm s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)
1775 !if (CS%id_hf_CAu > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_u,IsdB,IedB,jsd,jed,nz)
1776
1777 !CS%id_hf_CAv = register_diag_field('ocean_model', 'hf_CAv', diag%axesCvL, Time, &
1778 ! 'Fractional Thickness-weighted Meridional Coriolis and Advective Acceleration', &
1779 ! 'm s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)
1780 !if (CS%id_hf_CAv > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz)
1781
1782 cs%id_hf_PFu_2d = register_diag_field('ocean_model', 'hf_PFu_2d', diag%axesCu1, time, &
1783 'Depth-sum Fractional Thickness-weighted Zonal Pressure Force Acceleration', &
1784 'm s-2', conversion=us%L_T2_to_m_s2)
1785 if (cs%id_hf_PFu_2d > 0) call safe_alloc_ptr(cs%ADp%diag_hfrac_u,isdb,iedb,jsd,jed,nz)
1786
1787 cs%id_hf_PFv_2d = register_diag_field('ocean_model', 'hf_PFv_2d', diag%axesCv1, time, &
1788 'Depth-sum Fractional Thickness-weighted Meridional Pressure Force Acceleration', &
1789 'm s-2', conversion=us%L_T2_to_m_s2)
1790 if (cs%id_hf_PFv_2d > 0) call safe_alloc_ptr(cs%ADp%diag_hfrac_v,isd,ied,jsdb,jedb,nz)
1791
1792 cs%id_h_PFu = register_diag_field('ocean_model', 'h_PFu', diag%axesCuL, time, &
1793 'Thickness Multiplied Zonal Pressure Force Acceleration', &
1794 'm2 s-2', conversion=gv%H_to_m*us%L_T2_to_m_s2)
1795 if (cs%id_h_PFu > 0) call safe_alloc_ptr(cs%ADp%diag_hu,isdb,iedb,jsd,jed,nz)
1796
1797 cs%id_h_PFv = register_diag_field('ocean_model', 'h_PFv', diag%axesCvL, time, &
1798 'Thickness Multiplied Meridional Pressure Force Acceleration', &
1799 'm2 s-2', conversion=gv%H_to_m*us%L_T2_to_m_s2)
1800 if (cs%id_h_PFv > 0) call safe_alloc_ptr(cs%ADp%diag_hv,isd,ied,jsdb,jedb,nz)
1801
1802 cs%id_intz_PFu_2d = register_diag_field('ocean_model', 'intz_PFu_2d', diag%axesCu1, time, &
1803 'Depth-integral of Zonal Pressure Force Acceleration', &
1804 'm2 s-2', conversion=gv%H_to_m*us%L_T2_to_m_s2)
1805 if (cs%id_intz_PFu_2d > 0) call safe_alloc_ptr(cs%ADp%diag_hu,isdb,iedb,jsd,jed,nz)
1806
1807 cs%id_intz_PFv_2d = register_diag_field('ocean_model', 'intz_PFv_2d', diag%axesCv1, time, &
1808 'Depth-integral of Meridional Pressure Force Acceleration', &
1809 'm2 s-2', conversion=gv%H_to_m*us%L_T2_to_m_s2)
1810 if (cs%id_intz_PFv_2d > 0) call safe_alloc_ptr(cs%ADp%diag_hv,isd,ied,jsdb,jedb,nz)
1811
1812 cs%id_hf_CAu_2d = register_diag_field('ocean_model', 'hf_CAu_2d', diag%axesCu1, time, &
1813 'Depth-sum Fractional Thickness-weighted Zonal Coriolis and Advective Acceleration', &
1814 'm s-2', conversion=us%L_T2_to_m_s2)
1815 if (cs%id_hf_CAu_2d > 0) call safe_alloc_ptr(cs%ADp%diag_hfrac_u,isdb,iedb,jsd,jed,nz)
1816
1817 cs%id_hf_CAv_2d = register_diag_field('ocean_model', 'hf_CAv_2d', diag%axesCv1, time, &
1818 'Depth-sum Fractional Thickness-weighted Meridional Coriolis and Advective Acceleration', &
1819 'm s-2', conversion=us%L_T2_to_m_s2)
1820 if (cs%id_hf_CAv_2d > 0) call safe_alloc_ptr(cs%ADp%diag_hfrac_v,isd,ied,jsdb,jedb,nz)
1821
1822 cs%id_h_CAu = register_diag_field('ocean_model', 'h_CAu', diag%axesCuL, time, &
1823 'Thickness Multiplied Zonal Coriolis and Advective Acceleration', &
1824 'm2 s-2', conversion=gv%H_to_m*us%L_T2_to_m_s2)
1825 if (cs%id_h_CAu > 0) call safe_alloc_ptr(cs%ADp%diag_hu,isdb,iedb,jsd,jed,nz)
1826
1827 cs%id_h_CAv = register_diag_field('ocean_model', 'h_CAv', diag%axesCvL, time, &
1828 'Thickness Multiplied Meridional Coriolis and Advective Acceleration', &
1829 'm2 s-2', conversion=gv%H_to_m*us%L_T2_to_m_s2)
1830 if (cs%id_h_CAv > 0) call safe_alloc_ptr(cs%ADp%diag_hv,isd,ied,jsdb,jedb,nz)
1831
1832 cs%id_intz_CAu_2d = register_diag_field('ocean_model', 'intz_CAu_2d', diag%axesCu1, time, &
1833 'Depth-integral of Zonal Coriolis and Advective Acceleration', &
1834 'm2 s-2', conversion=gv%H_to_m*us%L_T2_to_m_s2)
1835 if (cs%id_intz_CAu_2d > 0) call safe_alloc_ptr(cs%ADp%diag_hu,isdb,iedb,jsd,jed,nz)
1836
1837 cs%id_intz_CAv_2d = register_diag_field('ocean_model', 'intz_CAv_2d', diag%axesCv1, time, &
1838 'Depth-integral of Meridional Coriolis and Advective Acceleration', &
1839 'm2 s-2', conversion=gv%H_to_m*us%L_T2_to_m_s2)
1840 if (cs%id_intz_CAv_2d > 0) call safe_alloc_ptr(cs%ADp%diag_hv,isd,ied,jsdb,jedb,nz)
1841
1842 cs%id_uav = register_diag_field('ocean_model', 'uav', diag%axesCuL, time, &
1843 'Barotropic-step Averaged Zonal Velocity', 'm s-1', conversion=us%L_T_to_m_s)
1844 cs%id_vav = register_diag_field('ocean_model', 'vav', diag%axesCvL, time, &
1845 'Barotropic-step Averaged Meridional Velocity', 'm s-1', conversion=us%L_T_to_m_s)
1846
1847 cs%id_u_BT_accel = register_diag_field('ocean_model', 'u_BT_accel', diag%axesCuL, time, &
1848 'Barotropic Anomaly Zonal Acceleration', 'm s-2', conversion=us%L_T2_to_m_s2)
1849 cs%id_v_BT_accel = register_diag_field('ocean_model', 'v_BT_accel', diag%axesCvL, time, &
1850 'Barotropic Anomaly Meridional Acceleration', 'm s-2', conversion=us%L_T2_to_m_s2)
1851
1852 !CS%id_hf_u_BT_accel = register_diag_field('ocean_model', 'hf_u_BT_accel', diag%axesCuL, Time, &
1853 ! 'Fractional Thickness-weighted Barotropic Anomaly Zonal Acceleration', &
1854 ! 'm s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)
1855 !if (CS%id_hf_u_BT_accel > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_u,IsdB,IedB,jsd,jed,nz)
1856
1857 !CS%id_hf_v_BT_accel = register_diag_field('ocean_model', 'hf_v_BT_accel', diag%axesCvL, Time, &
1858 ! 'Fractional Thickness-weighted Barotropic Anomaly Meridional Acceleration', &
1859 ! 'm s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)
1860 !if (CS%id_hf_v_BT_accel > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz)
1861
1862 cs%id_hf_u_BT_accel_2d = register_diag_field('ocean_model', 'hf_u_BT_accel_2d', diag%axesCu1, time, &
1863 'Depth-sum Fractional Thickness-weighted Barotropic Anomaly Zonal Acceleration', &
1864 'm s-2', conversion=us%L_T2_to_m_s2)
1865 if (cs%id_hf_u_BT_accel_2d > 0) call safe_alloc_ptr(cs%ADp%diag_hfrac_u,isdb,iedb,jsd,jed,nz)
1866
1867 cs%id_hf_v_BT_accel_2d = register_diag_field('ocean_model', 'hf_v_BT_accel_2d', diag%axesCv1, time, &
1868 'Depth-sum Fractional Thickness-weighted Barotropic Anomaly Meridional Acceleration', &
1869 'm s-2', conversion=us%L_T2_to_m_s2)
1870 if (cs%id_hf_v_BT_accel_2d > 0) call safe_alloc_ptr(cs%ADp%diag_hfrac_v,isd,ied,jsdb,jedb,nz)
1871
1872 cs%id_h_u_BT_accel = register_diag_field('ocean_model', 'h_u_BT_accel', diag%axesCuL, time, &
1873 'Thickness Multiplied Barotropic Anomaly Zonal Acceleration', &
1874 'm2 s-2', conversion=gv%H_to_m*us%L_T2_to_m_s2)
1875 if (cs%id_h_u_BT_accel > 0) call safe_alloc_ptr(cs%ADp%diag_hu,isdb,iedb,jsd,jed,nz)
1876
1877 cs%id_h_v_BT_accel = register_diag_field('ocean_model', 'h_v_BT_accel', diag%axesCvL, time, &
1878 'Thickness Multiplied Barotropic Anomaly Meridional Acceleration', &
1879 'm2 s-2', conversion=gv%H_to_m*us%L_T2_to_m_s2)
1880 if (cs%id_h_v_BT_accel > 0) call safe_alloc_ptr(cs%ADp%diag_hv,isd,ied,jsdb,jedb,nz)
1881
1882 cs%id_intz_u_BT_accel_2d = register_diag_field('ocean_model', 'intz_u_BT_accel_2d', diag%axesCu1, time, &
1883 'Depth-integral of Barotropic Anomaly Zonal Acceleration', &
1884 'm2 s-2', conversion=gv%H_to_m*us%L_T2_to_m_s2)
1885 if (cs%id_intz_u_BT_accel_2d > 0) call safe_alloc_ptr(cs%ADp%diag_hu,isdb,iedb,jsd,jed,nz)
1886
1887 cs%id_intz_v_BT_accel_2d = register_diag_field('ocean_model', 'intz_v_BT_accel_2d', diag%axesCv1, time, &
1888 'Depth-integral of Barotropic Anomaly Meridional Acceleration', &
1889 'm2 s-2', conversion=gv%H_to_m*us%L_T2_to_m_s2)
1890 if (cs%id_intz_v_BT_accel_2d > 0) call safe_alloc_ptr(cs%ADp%diag_hv,isd,ied,jsdb,jedb,nz)
1891
1892 cs%id_PFu_visc_rem = register_diag_field('ocean_model', 'PFu_visc_rem', diag%axesCuL, time, &
1893 'Zonal Pressure Force Acceleration multiplied by the viscous remnant', &
1894 'm s-2', conversion=us%L_T2_to_m_s2)
1895 if (cs%id_PFu_visc_rem > 0) call safe_alloc_ptr(cs%ADp%visc_rem_u,isdb,iedb,jsd,jed,nz)
1896 cs%id_PFv_visc_rem = register_diag_field('ocean_model', 'PFv_visc_rem', diag%axesCvL, time, &
1897 'Meridional Pressure Force Acceleration multiplied by the viscous remnant', &
1898 'm s-2', conversion=us%L_T2_to_m_s2)
1899 if (cs%id_PFv_visc_rem > 0) call safe_alloc_ptr(cs%ADp%visc_rem_v,isd,ied,jsdb,jedb,nz)
1900
1901 cs%id_CAu_visc_rem = register_diag_field('ocean_model', 'CAu_visc_rem', diag%axesCuL, time, &
1902 'Zonal Coriolis and Advective Acceleration multiplied by the viscous remnant', &
1903 'm s-2', conversion=us%L_T2_to_m_s2)
1904 if (cs%id_CAu_visc_rem > 0) call safe_alloc_ptr(cs%ADp%visc_rem_u,isdb,iedb,jsd,jed,nz)
1905 cs%id_CAv_visc_rem = register_diag_field('ocean_model', 'CAv_visc_rem', diag%axesCvL, time, &
1906 'Meridional Coriolis and Advective Acceleration multiplied by the viscous remnant', &
1907 'm s-2', conversion=us%L_T2_to_m_s2)
1908 if (cs%id_CAv_visc_rem > 0) call safe_alloc_ptr(cs%ADp%visc_rem_v,isd,ied,jsdb,jedb,nz)
1909
1910 cs%id_u_BT_accel_visc_rem = register_diag_field('ocean_model', 'u_BT_accel_visc_rem', diag%axesCuL, time, &
1911 'Barotropic Anomaly Zonal Acceleration multiplied by the viscous remnant', &
1912 'm s-2', conversion=us%L_T2_to_m_s2)
1913 if (cs%id_u_BT_accel_visc_rem > 0) call safe_alloc_ptr(cs%ADp%visc_rem_u,isdb,iedb,jsd,jed,nz)
1914 cs%id_v_BT_accel_visc_rem = register_diag_field('ocean_model', 'v_BT_accel_visc_rem', diag%axesCvL, time, &
1915 'Barotropic Anomaly Meridional Acceleration multiplied by the viscous remnant', &
1916 'm s-2', conversion=us%L_T2_to_m_s2)
1917 if (cs%id_v_BT_accel_visc_rem > 0) call safe_alloc_ptr(cs%ADp%visc_rem_v,isd,ied,jsdb,jedb,nz)
1918
1919 id_clock_cor = cpu_clock_id('(Ocean Coriolis & mom advection)', grain=clock_module)
1920 id_clock_continuity = cpu_clock_id('(Ocean continuity equation)', grain=clock_module)
1921 id_clock_pres = cpu_clock_id('(Ocean pressure force)', grain=clock_module)
1922 id_clock_vertvisc = cpu_clock_id('(Ocean vertical viscosity)', grain=clock_module)
1923 id_clock_horvisc = cpu_clock_id('(Ocean horizontal viscosity)', grain=clock_module)
1924 id_clock_mom_update = cpu_clock_id('(Ocean momentum increments)', grain=clock_module)
1925 id_clock_pass = cpu_clock_id('(Ocean message passing)', grain=clock_module)
1926 id_clock_btcalc = cpu_clock_id('(Ocean barotropic mode calc)', grain=clock_module)
1927 id_clock_btstep = cpu_clock_id('(Ocean barotropic mode stepping)', grain=clock_module)
1928 id_clock_btforce = cpu_clock_id('(Ocean barotropic forcing calc)', grain=clock_module)
1929
1930end subroutine initialize_dyn_split_rk2
1931
1932
1933!> Close the dyn_split_RK2 module
1934subroutine end_dyn_split_rk2(CS)
1935 type(mom_dyn_split_rk2_cs), pointer :: cs !< module control structure
1936
1937 call barotropic_end(cs%barotropic_CSp)
1938
1939 call vertvisc_end(cs%vertvisc_CSp)
1940 deallocate(cs%vertvisc_CSp)
1941
1942 call hor_visc_end(cs%hor_visc)
1943 if (cs%calculate_SAL) call sal_end(cs%SAL_CSp)
1944 if (cs%use_tides) call tidal_forcing_end(cs%tides_CSp)
1945 call coriolisadv_end(cs%CoriolisAdv)
1946
1947 dealloc_(cs%diffu) ; dealloc_(cs%diffv)
1948 dealloc_(cs%CAu) ; dealloc_(cs%CAv)
1949 dealloc_(cs%CAu_pred) ; dealloc_(cs%CAv_pred)
1950 dealloc_(cs%PFu) ; dealloc_(cs%PFv)
1951
1952 if (associated(cs%taux_bot)) deallocate(cs%taux_bot)
1953 if (associated(cs%tauy_bot)) deallocate(cs%tauy_bot)
1954
1955 dealloc_(cs%uhbt) ; dealloc_(cs%vhbt)
1956 dealloc_(cs%u_accel_bt) ; dealloc_(cs%v_accel_bt)
1957 dealloc_(cs%visc_rem_u) ; dealloc_(cs%visc_rem_v)
1958
1959 dealloc_(cs%eta) ; dealloc_(cs%eta_PF) ; dealloc_(cs%pbce)
1960 dealloc_(cs%h_av) ; dealloc_(cs%u_av) ; dealloc_(cs%v_av)
1961
1962 call dealloc_bt_cont_type(cs%BT_cont)
1963 deallocate(cs%AD_pred)
1964
1965 deallocate(cs)
1966end subroutine end_dyn_split_rk2
1967
1968
1969!> \namespace mom_dynamics_split_rk2
1970!!
1971!! This file time steps the adiabatic dynamic core by splitting
1972!! between baroclinic and barotropic modes. It uses a pseudo-second order
1973!! Runge-Kutta time stepping scheme for the baroclinic momentum
1974!! equation and a forward-backward coupling between the baroclinic
1975!! momentum and continuity equations. This split time-stepping
1976!! scheme is described in detail in Hallberg (JCP, 1997). Additional
1977!! issues related to exact tracer conservation and how to
1978!! ensure consistency between the barotropic and layered estimates
1979!! of the free surface height are described in Hallberg and
1980!! Adcroft (Ocean Modelling, 2009). This was the time stepping code
1981!! that is used for most GOLD applications, including GFDL's ESM2G
1982!! Earth system model, and all of the examples provided with the
1983!! MOM code (although several of these solutions are routinely
1984!! verified by comparison with the slower unsplit schemes).
1985!!
1986!! The subroutine step_MOM_dyn_split_RK2 actually does the time
1987!! stepping, while register_restarts_dyn_split_RK2 sets the fields
1988!! that are found in a full restart file with this scheme, and
1989!! initialize_dyn_split_RK2 initializes the cpu clocks that are
1990!! used in this module. For largely historical reasons, this module
1991!! does not have its own control structure, but shares the same
1992!! control structure with MOM.F90 and the other MOM_dynamics_...
1993!! modules.
1994
1995end module mom_dynamics_split_rk2