MOM6
MOM_dynamics_split_RK2.F90
1 !> Time step the adiabatic dynamic core of MOM using RK2 method.
2 module mom_dynamics_split_rk2
3 
4 ! This file is part of MOM6. See LICENSE.md for the license.
5 
6 use mom_variables, only : vertvisc_type, thermo_var_ptrs
7 use mom_variables, only : bt_cont_type, alloc_bt_cont_type, dealloc_bt_cont_type
8 use mom_variables, only : accel_diag_ptrs, ocean_internal_state, cont_diag_ptrs
9 use mom_forcing_type, only : mech_forcing
10 
11 use mom_checksum_packages, only : mom_thermo_chksum, mom_state_chksum, mom_accel_chksum
12 use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
13 use mom_cpu_clock, only : clock_component, clock_subcomponent
14 use mom_cpu_clock, only : clock_module_driver, clock_module, clock_routine
15 use mom_diag_mediator, only : diag_mediator_init, enable_averages
16 use mom_diag_mediator, only : disable_averaging, post_data, safe_alloc_ptr
17 use mom_diag_mediator, only : register_diag_field, register_static_field
18 use mom_diag_mediator, only : set_diag_mediator_grid, diag_ctrl, diag_update_remap_grids
19 use mom_domains, only : mom_domains_init
20 use mom_domains, only : to_south, to_west, to_all, cgrid_ne, scalar_pair
21 use mom_domains, only : to_north, to_east, omit_corners
22 use mom_domains, only : create_group_pass, do_group_pass, group_pass_type
23 use mom_domains, only : start_group_pass, complete_group_pass, pass_var
24 use mom_debugging, only : hchksum, uvchksum
25 use mom_error_handler, only : mom_error, mom_mesg, fatal, warning, is_root_pe
26 use mom_error_handler, only : mom_set_verbosity, calltree_showquery
27 use mom_error_handler, only : calltree_enter, calltree_leave, calltree_waypoint
28 use mom_file_parser, only : get_param, log_version, param_file_type
29 use mom_get_input, only : directories
30 use mom_io, only : mom_io_init, vardesc, var_desc
31 use mom_restart, only : register_restart_field, register_restart_pair
32 use mom_restart, only : query_initialized, save_restart
33 use mom_restart, only : restart_init, is_new_run, mom_restart_cs
34 use mom_time_manager, only : time_type, time_type_to_real, operator(+)
35 use mom_time_manager, only : operator(-), operator(>), operator(*), operator(/)
36 
37 use mom_ale, only : ale_cs
38 use mom_barotropic, only : barotropic_init, btstep, btcalc, bt_mass_source
39 use mom_barotropic, only : register_barotropic_restarts, set_dtbt, barotropic_cs
40 use mom_boundary_update, only : update_obc_data, update_obc_cs
41 use mom_continuity, only : continuity, continuity_init, continuity_cs
42 use mom_continuity, only : continuity_stencil
43 use mom_coriolisadv, only : coradcalc, coriolisadv_init, coriolisadv_cs
44 use mom_debugging, only : check_redundant
45 use mom_grid, only : ocean_grid_type
46 use mom_hor_index, only : hor_index_type
47 use mom_hor_visc, only : horizontal_viscosity, hor_visc_init, hor_visc_cs
48 use mom_interface_heights, only : find_eta
49 use mom_lateral_mixing_coeffs, only : varmix_cs
50 use mom_meke_types, only : meke_type
51 use mom_open_boundary, only : ocean_obc_type, radiation_open_bdry_conds
52 use mom_open_boundary, only : open_boundary_zero_normal_flow
53 use mom_open_boundary, only : open_boundary_test_extern_h, update_obc_ramp
54 use mom_pressureforce, only : pressureforce, pressureforce_init, pressureforce_cs
55 use mom_set_visc, only : set_viscous_ml, set_visc_cs
56 use mom_thickness_diffuse, only : thickness_diffuse_cs
57 use mom_tidal_forcing, only : tidal_forcing_init, tidal_forcing_cs
58 use mom_unit_scaling, only : unit_scale_type
59 use mom_vert_friction, only : vertvisc, vertvisc_coef, vertvisc_remnant
60 use mom_vert_friction, only : vertvisc_init, vertvisc_cs
61 use mom_vert_friction, only : updatecfltruncationvalue
62 use mom_verticalgrid, only : verticalgrid_type, get_thickness_units
63 use mom_verticalgrid, only : get_flux_units, get_tr_flux_units
64 use mom_wave_interface, only: wave_parameters_cs
65 
66 implicit none ; private
67 
68 #include <MOM_memory.h>
69 
70 !> MOM_dynamics_split_RK2 module control structure
71 type, public :: mom_dyn_split_rk2_cs ; private
72  real allocable_, dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: &
73  cau, & !< CAu = f*v - u.grad(u) [L T-2 ~> m s-2]
74  pfu, & !< PFu = -dM/dx [L T-2 ~> m s-2]
75  diffu !< Zonal acceleration due to convergence of the along-isopycnal stress tensor [L T-2 ~> m s-2]
76 
77  real allocable_, dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: &
78  cav, & !< CAv = -f*u - u.grad(v) [L T-2 ~> m s-2]
79  pfv, & !< PFv = -dM/dy [L T-2 ~> m s-2]
80  diffv !< Meridional acceleration due to convergence of the along-isopycnal stress tensor [L T-2 ~> m s-2]
81 
82  real allocable_, dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: visc_rem_u
83  !< Both the fraction of the zonal momentum originally in a
84  !! layer that remains after a time-step of viscosity, and the
85  !! fraction of a time-step worth of a barotropic acceleration
86  !! that a layer experiences after viscosity is applied.
87  !! Nondimensional between 0 (at the bottom) and 1 (far above).
88  real allocable_, dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: u_accel_bt
89  !< The zonal layer accelerations due to the difference between
90  !! the barotropic accelerations and the baroclinic accelerations
91  !! that were fed into the barotopic calculation [L T-2 ~> m s-2]
92  real allocable_, dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: visc_rem_v
93  !< Both the fraction of the meridional momentum originally in
94  !! a layer that remains after a time-step of viscosity, and the
95  !! fraction of a time-step worth of a barotropic acceleration
96  !! that a layer experiences after viscosity is applied.
97  !! Nondimensional between 0 (at the bottom) and 1 (far above).
98  real allocable_, dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: v_accel_bt
99  !< The meridional layer accelerations due to the difference between
100  !! the barotropic accelerations and the baroclinic accelerations
101  !! that were fed into the barotopic calculation [L T-2 ~> m s-2]
102 
103  ! The following variables are only used with the split time stepping scheme.
104  real allocable_, dimension(NIMEM_,NJMEM_) :: eta !< Instantaneous free surface height (in Boussinesq
105  !! mode) or column mass anomaly (in non-Boussinesq
106  !! mode) [H ~> m or kg m-2]
107  real allocable_, dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: u_av !< layer x-velocity with vertical mean replaced by
108  !! time-mean barotropic velocity over a baroclinic
109  !! timestep [L T-1 ~> m s-1]
110  real allocable_, dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: v_av !< layer y-velocity with vertical mean replaced by
111  !! time-mean barotropic velocity over a baroclinic
112  !! timestep [L T-1 ~> m s-1]
113  real allocable_, dimension(NIMEM_,NJMEM_,NKMEM_) :: h_av !< arithmetic mean of two successive layer
114  !! thicknesses [H ~> m or kg m-2]
115  real allocable_, dimension(NIMEM_,NJMEM_) :: eta_pf !< instantaneous SSH used in calculating PFu and
116  !! PFv [H ~> m or kg m-2]
117  real allocable_, dimension(NIMEMB_PTR_,NJMEM_) :: uhbt !< average x-volume or mass flux determined by the
118  !! barotropic solver [H L2 T-1 ~> m3 s-1 or kg s-1].
119  !! uhbt is roughly equal to the vertical sum of uh.
120  real allocable_, dimension(NIMEM_,NJMEMB_PTR_) :: vhbt !< average y-volume or mass flux determined by the
121  !! barotropic solver [H L2 T-1 ~> m3 s-1 or kg s-1].
122  !! vhbt is roughly equal to vertical sum of vh.
123  real allocable_, dimension(NIMEM_,NJMEM_,NKMEM_) :: pbce !< pbce times eta gives the baroclinic pressure
124  !! anomaly in each layer due to free surface height
125  !! anomalies [L2 H-1 T-2 ~> m s-2 or m4 kg-1 s-2].
126 
127  real, pointer, dimension(:,:) :: taux_bot => null() !< frictional x-bottom stress from the ocean
128  !! to the seafloor [R L Z T-2 ~> Pa]
129  real, pointer, dimension(:,:) :: tauy_bot => null() !< frictional y-bottom stress from the ocean
130  !! to the seafloor [R L Z T-2 ~> Pa]
131  type(bt_cont_type), pointer :: bt_cont => null() !< A structure with elements that describe the
132  !! effective summed open face areas as a function
133  !! of barotropic flow.
134 
135  ! This is to allow the previous, velocity-based coupling with between the
136  ! baroclinic and barotropic modes.
137  logical :: bt_use_layer_fluxes !< If true, use the summed layered fluxes plus
138  !! an adjustment due to a changed barotropic
139  !! velocity in the barotropic continuity equation.
140  logical :: split_bottom_stress !< If true, provide the bottom stress
141  !! calculated by the vertical viscosity to the
142  !! barotropic solver.
143  logical :: calc_dtbt !< If true, calculate the barotropic time-step
144  !! dynamically.
145 
146  real :: be !< A nondimensional number from 0.5 to 1 that controls
147  !! the backward weighting of the time stepping scheme.
148  real :: begw !< A nondimensional number from 0 to 1 that controls
149  !! the extent to which the treatment of gravity waves
150  !! is forward-backward (0) or simulated backward
151  !! Euler (1). 0 is almost always used.
152  logical :: debug !< If true, write verbose checksums for debugging purposes.
153  logical :: debug_obc !< If true, do debugging calls for open boundary conditions.
154 
155  logical :: module_is_initialized = .false. !< Record whether this mouled has been initialzed.
156 
157  !>@{ Diagnostic IDs
158  integer :: id_uh = -1, id_vh = -1
159  integer :: id_umo = -1, id_vmo = -1
160  integer :: id_umo_2d = -1, id_vmo_2d = -1
161  integer :: id_pfu = -1, id_pfv = -1
162  integer :: id_cau = -1, id_cav = -1
163  ! integer :: id_hf_PFu = -1, id_hf_PFv = -1
164  integer :: id_hf_pfu_2d = -1, id_hf_pfv_2d = -1
165  ! integer :: id_hf_CAu = -1, id_hf_CAv = -1
166  integer :: id_hf_cau_2d = -1, id_hf_cav_2d = -1
167 
168  ! Split scheme only.
169  integer :: id_uav = -1, id_vav = -1
170  integer :: id_u_bt_accel = -1, id_v_bt_accel = -1
171  ! integer :: id_hf_u_BT_accel = -1, id_hf_v_BT_accel = -1
172  integer :: id_hf_u_bt_accel_2d = -1, id_hf_v_bt_accel_2d = -1
173  !>@}
174 
175  type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the
176  !! timing of diagnostic output.
177  type(accel_diag_ptrs), pointer :: adp !< A structure pointing to the various
178  !! accelerations in the momentum equations,
179  !! which can later be used to calculate
180  !! derived diagnostics like energy budgets.
181  type(cont_diag_ptrs), pointer :: cdp !< A structure with pointers to various
182  !! terms in the continuity equations,
183  !! which can later be used to calculate
184  !! derived diagnostics like energy budgets.
185 
186  ! The remainder of the structure points to child subroutines' control structures.
187  !> A pointer to the horizontal viscosity control structure
188  type(hor_visc_cs), pointer :: hor_visc_csp => null()
189  !> A pointer to the continuity control structure
190  type(continuity_cs), pointer :: continuity_csp => null()
191  !> A pointer to the CoriolisAdv control structure
192  type(coriolisadv_cs), pointer :: coriolisadv_csp => null()
193  !> A pointer to the PressureForce control structure
194  type(pressureforce_cs), pointer :: pressureforce_csp => null()
195  !> A pointer to the barotropic stepping control structure
196  type(barotropic_cs), pointer :: barotropic_csp => null()
197  !> A pointer to a structure containing interface height diffusivities
198  type(thickness_diffuse_cs), pointer :: thickness_diffuse_csp => null()
199  !> A pointer to the vertical viscosity control structure
200  type(vertvisc_cs), pointer :: vertvisc_csp => null()
201  !> A pointer to the set_visc control structure
202  type(set_visc_cs), pointer :: set_visc_csp => null()
203  !> A pointer to the tidal forcing control structure
204  type(tidal_forcing_cs), pointer :: tides_csp => null()
205  !> A pointer to the ALE control structure.
206  type(ale_cs), pointer :: ale_csp => null()
207 
208  type(ocean_obc_type), pointer :: obc => null() !< A pointer to an open boundary
209  !! condition type that specifies whether, where, and what open boundary
210  !! conditions are used. If no open BCs are used, this pointer stays
211  !! nullified. Flather OBCs use open boundary_CS as well.
212  !> A pointer to the update_OBC control structure
213  type(update_obc_cs), pointer :: update_obc_csp => null()
214 
215  type(group_pass_type) :: pass_eta !< Structure for group halo pass
216  type(group_pass_type) :: pass_visc_rem !< Structure for group halo pass
217  type(group_pass_type) :: pass_uvp !< Structure for group halo pass
218  type(group_pass_type) :: pass_hp_uv !< Structure for group halo pass
219  type(group_pass_type) :: pass_uv !< Structure for group halo pass
220  type(group_pass_type) :: pass_h !< Structure for group halo pass
221  type(group_pass_type) :: pass_av_uvh !< Structure for group halo pass
222 
223 end type mom_dyn_split_rk2_cs
224 
225 
226 public step_mom_dyn_split_rk2
227 public register_restarts_dyn_split_rk2
228 public initialize_dyn_split_rk2
229 public end_dyn_split_rk2
230 
231 !>@{ CPU time clock IDs
232 integer :: id_clock_cor, id_clock_pres, id_clock_vertvisc
233 integer :: id_clock_horvisc, id_clock_mom_update
234 integer :: id_clock_continuity, id_clock_thick_diff
235 integer :: id_clock_btstep, id_clock_btcalc, id_clock_btforce
236 integer :: id_clock_pass, id_clock_pass_init
237 !>@}
238 
239 contains
240 
241 !> RK2 splitting for time stepping MOM adiabatic dynamics
242 subroutine step_mom_dyn_split_rk2(u, v, h, tv, visc, Time_local, dt, forces, p_surf_begin, p_surf_end, &
243  uh, vh, uhtr, vhtr, eta_av, G, GV, US, CS, calc_dtbt, VarMix, &
244  MEKE, thickness_diffuse_CSp, Waves)
245  type(ocean_grid_type), intent(inout) :: g !< ocean grid structure
246  type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
247  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
248  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
249  target, intent(inout) :: u !< zonal velocity [L T-1 ~> m s-1]
250  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
251  target, intent(inout) :: v !< merid velocity [L T-1 ~> m s-1]
252  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
253  intent(inout) :: h !< layer thickness [H ~> m or kg m-2]
254  type(thermo_var_ptrs), intent(in) :: tv !< thermodynamic type
255  type(vertvisc_type), intent(inout) :: visc !< vertical visc, bottom drag, and related
256  type(time_type), intent(in) :: time_local !< model time at end of time step
257  real, intent(in) :: dt !< time step [T ~> s]
258  type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces
259  real, dimension(:,:), pointer :: p_surf_begin !< surf pressure at the start of this dynamic
260  !! time step [R L2 T-2 ~> Pa]
261  real, dimension(:,:), pointer :: p_surf_end !< surf pressure at the end of this dynamic
262  !! time step [R L2 T-2 ~> Pa]
263  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
264  target, intent(inout) :: uh !< zonal volume/mass transport
265  !! [H L2 T-1 ~> m3 s-1 or kg s-1]
266  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
267  target, intent(inout) :: vh !< merid volume/mass transport
268  !! [H L2 T-1 ~> m3 s-1 or kg s-1]
269  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
270  intent(inout) :: uhtr !< accumulatated zonal volume/mass transport
271  !! since last tracer advection [H L2 ~> m3 or kg]
272  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
273  intent(inout) :: vhtr !< accumulatated merid volume/mass transport
274  !! since last tracer advection [H L2 ~> m3 or kg]
275  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: eta_av !< free surface height or column mass time
276  !! averaged over time step [H ~> m or kg m-2]
277  type(mom_dyn_split_rk2_cs), pointer :: cs !< module control structure
278  logical, intent(in) :: calc_dtbt !< if true, recalculate barotropic time step
279  type(varmix_cs), pointer :: varmix !< specify the spatially varying viscosities
280  type(meke_type), pointer :: meke !< related to mesoscale eddy kinetic energy param
281  type(thickness_diffuse_cs), pointer :: thickness_diffuse_csp !< Pointer to a structure containing
282  !! interface height diffusivities
283  type(wave_parameters_cs), optional, pointer :: waves !< A pointer to a structure containing
284  !! fields related to the surface wave conditions
285 
286  ! local variables
287  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: up ! Predicted zonal velocity [L T-1 ~> m s-1].
288  real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: vp ! Predicted meridional velocity [L T-1 ~> m s-1].
289  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: hp ! Predicted thickness [H ~> m or kg m-2].
290 
291  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: u_bc_accel
292  real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: v_bc_accel
293  ! u_bc_accel and v_bc_accel are the summed baroclinic accelerations of each
294  ! layer calculated by the non-barotropic part of the model [L T-2 ~> m s-2].
295 
296  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), target :: uh_in
297  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), target :: vh_in
298  ! uh_in and vh_in are the zonal or meridional mass transports that would be
299  ! obtained using the initial velocities [H L2 T-1 ~> m3 s-1 or kg s-1].
300 
301  real, dimension(SZIB_(G),SZJ_(G)) :: uhbt_out
302  real, dimension(SZI_(G),SZJB_(G)) :: vhbt_out
303  ! uhbt_out and vhbt_out are the vertically summed transports from the
304  ! barotropic solver based on its final velocities [H m2 s-1 ~> m3 s-1 or kg s-1].
305 
306  real, dimension(SZI_(G),SZJ_(G)) :: eta_pred
307  ! eta_pred is the predictor value of the free surface height or column mass,
308  ! [H ~> m or kg m-2].
309 
310  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: u_old_rad_obc
311  real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: v_old_rad_obc
312  ! u_old_rad_OBC and v_old_rad_OBC are the starting velocities, which are
313  ! saved for use in the Flather open boundary condition code [L T-1 ~> m s-1].
314 
315  real :: pres_to_eta ! A factor that converts pressures to the units of eta
316  ! [H T2 R-1 L-2 ~> m Pa-1 or kg m-2 Pa-1]
317  real, pointer, dimension(:,:) :: &
318  p_surf => null(), eta_pf_start => null(), &
319  taux_bot => null(), tauy_bot => null(), &
320  eta => null()
321 
322  real, pointer, dimension(:,:,:) :: &
323  uh_ptr => null(), u_ptr => null(), vh_ptr => null(), v_ptr => null(), &
324  u_av, & ! The zonal velocity time-averaged over a time step [L T-1 ~> m s-1].
325  v_av, & ! The meridional velocity time-averaged over a time step [L T-1 ~> m s-1].
326  h_av ! The layer thickness time-averaged over a time step [H ~> m or kg m-2].
327 
328  ! real, allocatable, dimension(:,:,:) :: &
329  ! hf_PFu, hf_PFv, & ! Pressure force accel. x fract. thickness [L T-2 ~> m s-2].
330  ! hf_CAu, hf_CAv, & ! Coriolis force accel. x fract. thickness [L T-2 ~> m s-2].
331  ! hf_u_BT_accel, hf_v_BT_accel ! barotropic correction accel. x fract. thickness [L T-2 ~> m s-2].
332  ! 3D diagnostics hf_PFu etc. are commented because there is no clarity on proper remapping grid option.
333  ! The code is retained for degugging purposes in the future.
334 
335  real, allocatable, dimension(:,:) :: &
336  hf_pfu_2d, hf_pfv_2d, & ! Depth integeral of hf_PFu, hf_PFv [L T-2 ~> m s-2].
337  hf_cau_2d, hf_cav_2d, & ! Depth integeral of hf_CAu, hf_CAv [L T-2 ~> m s-2].
338  hf_u_bt_accel_2d, hf_v_bt_accel_2d ! Depth integeral of hf_u_BT_accel, hf_v_BT_accel
339 
340  real :: dt_pred ! The time step for the predictor part of the baroclinic time stepping [T ~> s].
341 
342  logical :: dyn_p_surf
343  logical :: bt_cont_bt_thick ! If true, use the BT_cont_type to estimate the
344  ! relative weightings of the layers in calculating
345  ! the barotropic accelerations.
346  !---For group halo pass
347  logical :: showcalltree, sym
348 
349  integer :: i, j, k, is, ie, js, je, isq, ieq, jsq, jeq, nz
350  integer :: cont_stencil
351  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
352  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
353  u_av => cs%u_av ; v_av => cs%v_av ; h_av => cs%h_av ; eta => cs%eta
354 
355  sym=.false.;if (g%Domain%symmetric) sym=.true. ! switch to include symmetric domain in checksums
356 
357  showcalltree = calltree_showquery()
358  if (showcalltree) call calltree_enter("step_MOM_dyn_split_RK2(), MOM_dynamics_split_RK2.F90")
359 
360  !$OMP parallel do default(shared)
361  do k=1,nz
362  do j=g%jsd,g%jed ; do i=g%isdB,g%iedB ; up(i,j,k) = 0.0 ; enddo ; enddo
363  do j=g%jsdB,g%jedB ; do i=g%isd,g%ied ; vp(i,j,k) = 0.0 ; enddo ; enddo
364  do j=g%jsd,g%jed ; do i=g%isd,g%ied ; hp(i,j,k) = h(i,j,k) ; enddo ; enddo
365  enddo
366 
367  ! Update CFL truncation value as function of time
368  call updatecfltruncationvalue(time_local, cs%vertvisc_CSp)
369 
370  if (cs%debug) then
371  call mom_state_chksum("Start predictor ", u, v, h, uh, vh, g, gv, us, symmetric=sym)
372  call check_redundant("Start predictor u ", u, v, g)
373  call check_redundant("Start predictor uh ", uh, vh, g)
374  endif
375 
376  dyn_p_surf = associated(p_surf_begin) .and. associated(p_surf_end)
377  if (dyn_p_surf) then
378  p_surf => p_surf_end
379  call safe_alloc_ptr(eta_pf_start,g%isd,g%ied,g%jsd,g%jed)
380  eta_pf_start(:,:) = 0.0
381  else
382  p_surf => forces%p_surf
383  endif
384 
385  if (associated(cs%OBC)) then
386  if (cs%debug_OBC) call open_boundary_test_extern_h(g, gv, cs%OBC, h)
387 
388  ! Update OBC ramp value as function of time
389  call update_obc_ramp(time_local, cs%OBC)
390 
391  do k=1,nz ; do j=g%jsd,g%jed ; do i=g%IsdB,g%IedB
392  u_old_rad_obc(i,j,k) = u_av(i,j,k)
393  enddo ; enddo ; enddo
394  do k=1,nz ; do j=g%JsdB,g%JedB ; do i=g%isd,g%ied
395  v_old_rad_obc(i,j,k) = v_av(i,j,k)
396  enddo ; enddo ; enddo
397  endif
398 
399  bt_cont_bt_thick = .false.
400  if (associated(cs%BT_cont)) bt_cont_bt_thick = &
401  (allocated(cs%BT_cont%h_u) .and. allocated(cs%BT_cont%h_v))
402 
403  if (cs%split_bottom_stress) then
404  taux_bot => cs%taux_bot ; tauy_bot => cs%tauy_bot
405  endif
406 
407  !--- begin set up for group halo pass
408 
409  cont_stencil = continuity_stencil(cs%continuity_CSp)
410  call cpu_clock_begin(id_clock_pass)
411  call create_group_pass(cs%pass_eta, eta, g%Domain, halo=1)
412  call create_group_pass(cs%pass_visc_rem, cs%visc_rem_u, cs%visc_rem_v, g%Domain, &
413  to_all+scalar_pair, cgrid_ne, halo=max(1,cont_stencil))
414  call create_group_pass(cs%pass_uvp, up, vp, g%Domain, halo=max(1,cont_stencil))
415  call create_group_pass(cs%pass_hp_uv, hp, g%Domain, halo=2)
416  call create_group_pass(cs%pass_hp_uv, u_av, v_av, g%Domain, halo=2)
417  call create_group_pass(cs%pass_hp_uv, uh(:,:,:), vh(:,:,:), g%Domain, halo=2)
418 
419  call create_group_pass(cs%pass_uv, u, v, g%Domain, halo=max(2,cont_stencil))
420  call create_group_pass(cs%pass_h, h, g%Domain, halo=max(2,cont_stencil))
421  call create_group_pass(cs%pass_av_uvh, u_av, v_av, g%Domain, halo=2)
422  call create_group_pass(cs%pass_av_uvh, uh(:,:,:), vh(:,:,:), g%Domain, halo=2)
423  call cpu_clock_end(id_clock_pass)
424  !--- end set up for group halo pass
425 
426 
427 ! PFu = d/dx M(h,T,S)
428 ! pbce = dM/deta
429  if (cs%begw == 0.0) call enable_averages(dt, time_local, cs%diag)
430  call cpu_clock_begin(id_clock_pres)
431  call pressureforce(h, tv, cs%PFu, cs%PFv, g, gv, us, cs%PressureForce_CSp, &
432  cs%ALE_CSp, p_surf, cs%pbce, cs%eta_PF)
433  if (dyn_p_surf) then
434  pres_to_eta = 1.0 / (gv%g_Earth * gv%H_to_RZ)
435  !$OMP parallel do default(shared)
436  do j=jsq,jeq+1 ; do i=isq,ieq+1
437  eta_pf_start(i,j) = cs%eta_PF(i,j) - pres_to_eta * (p_surf_begin(i,j) - p_surf_end(i,j))
438  enddo ; enddo
439  endif
440  call cpu_clock_end(id_clock_pres)
441  call disable_averaging(cs%diag)
442  if (showcalltree) call calltree_waypoint("done with PressureForce (step_MOM_dyn_split_RK2)")
443 
444  if (associated(cs%OBC)) then; if (cs%OBC%update_OBC) then
445  call update_obc_data(cs%OBC, g, gv, us, tv, h, cs%update_OBC_CSp, time_local)
446  endif; endif
447  if (associated(cs%OBC) .and. cs%debug_OBC) &
448  call open_boundary_zero_normal_flow(cs%OBC, g, cs%PFu, cs%PFv)
449 
450  if (g%nonblocking_updates) &
451  call start_group_pass(cs%pass_eta, g%Domain, clock=id_clock_pass)
452 
453 ! CAu = -(f+zeta_av)/h_av vh + d/dx KE_av
454  call cpu_clock_begin(id_clock_cor)
455  call coradcalc(u_av, v_av, h_av, uh, vh, cs%CAu, cs%CAv, cs%OBC, cs%ADp, &
456  g, gv, us, cs%CoriolisAdv_CSp)
457  call cpu_clock_end(id_clock_cor)
458  if (showcalltree) call calltree_waypoint("done with CorAdCalc (step_MOM_dyn_split_RK2)")
459 
460 ! u_bc_accel = CAu + PFu + diffu(u[n-1])
461  call cpu_clock_begin(id_clock_btforce)
462  !$OMP parallel do default(shared)
463  do k=1,nz
464  do j=js,je ; do i=isq,ieq
465  u_bc_accel(i,j,k) = (cs%CAu(i,j,k) + cs%PFu(i,j,k)) + cs%diffu(i,j,k)
466  enddo ; enddo
467  do j=jsq,jeq ; do i=is,ie
468  v_bc_accel(i,j,k) = (cs%CAv(i,j,k) + cs%PFv(i,j,k)) + cs%diffv(i,j,k)
469  enddo ; enddo
470  enddo
471  if (associated(cs%OBC)) then
472  call open_boundary_zero_normal_flow(cs%OBC, g, u_bc_accel, v_bc_accel)
473  endif
474  call cpu_clock_end(id_clock_btforce)
475 
476  if (cs%debug) then
477  call mom_accel_chksum("pre-btstep accel", cs%CAu, cs%CAv, cs%PFu, cs%PFv, &
478  cs%diffu, cs%diffv, g, gv, us, cs%pbce, u_bc_accel, v_bc_accel, &
479  symmetric=sym)
480  call check_redundant("pre-btstep CS%Ca ", cs%Cau, cs%Cav, g)
481  call check_redundant("pre-btstep CS%PF ", cs%PFu, cs%PFv, g)
482  call check_redundant("pre-btstep CS%diff ", cs%diffu, cs%diffv, g)
483  call check_redundant("pre-btstep u_bc_accel ", u_bc_accel, v_bc_accel, g)
484  endif
485 
486  call cpu_clock_begin(id_clock_vertvisc)
487  !$OMP parallel do default(shared)
488  do k=1,nz
489  do j=js,je ; do i=isq,ieq
490  up(i,j,k) = g%mask2dCu(i,j) * (u(i,j,k) + dt * u_bc_accel(i,j,k))
491  enddo ; enddo
492  do j=jsq,jeq ; do i=is,ie
493  vp(i,j,k) = g%mask2dCv(i,j) * (v(i,j,k) + dt * v_bc_accel(i,j,k))
494  enddo ; enddo
495  enddo
496 
497  call enable_averages(dt, time_local, cs%diag)
498  call set_viscous_ml(u, v, h, tv, forces, visc, dt, g, gv, us, &
499  cs%set_visc_CSp)
500  call disable_averaging(cs%diag)
501 
502  if (cs%debug) then
503  call uvchksum("before vertvisc: up", up, vp, g%HI, haloshift=0, symmetric=sym, scale=us%L_T_to_m_s)
504  endif
505  call vertvisc_coef(up, vp, h, forces, visc, dt, g, gv, us, cs%vertvisc_CSp, cs%OBC)
506  call vertvisc_remnant(visc, cs%visc_rem_u, cs%visc_rem_v, dt, g, gv, us, cs%vertvisc_CSp)
507  call cpu_clock_end(id_clock_vertvisc)
508  if (showcalltree) call calltree_waypoint("done with vertvisc_coef (step_MOM_dyn_split_RK2)")
509 
510 
511  call cpu_clock_begin(id_clock_pass)
512  if (g%nonblocking_updates) then
513  call complete_group_pass(cs%pass_eta, g%Domain)
514  call start_group_pass(cs%pass_visc_rem, g%Domain)
515  else
516  call do_group_pass(cs%pass_eta, g%Domain)
517  call do_group_pass(cs%pass_visc_rem, g%Domain)
518  endif
519  call cpu_clock_end(id_clock_pass)
520 
521  call cpu_clock_begin(id_clock_btcalc)
522  ! Calculate the relative layer weights for determining barotropic quantities.
523  if (.not.bt_cont_bt_thick) &
524  call btcalc(h, g, gv, cs%barotropic_CSp, obc=cs%OBC)
525  call bt_mass_source(h, eta, .true., g, gv, cs%barotropic_CSp)
526  call cpu_clock_end(id_clock_btcalc)
527 
528  if (g%nonblocking_updates) &
529  call complete_group_pass(cs%pass_visc_rem, g%Domain, clock=id_clock_pass)
530 
531 ! u_accel_bt = layer accelerations due to barotropic solver
532  if (associated(cs%BT_cont) .or. cs%BT_use_layer_fluxes) then
533  call cpu_clock_begin(id_clock_continuity)
534  call continuity(u, v, h, hp, uh_in, vh_in, dt, g, gv, us, cs%continuity_CSp, &
535  obc=cs%OBC, visc_rem_u=cs%visc_rem_u, visc_rem_v=cs%visc_rem_v, bt_cont=cs%BT_cont)
536  call cpu_clock_end(id_clock_continuity)
537  if (bt_cont_bt_thick) then
538  call btcalc(h, g, gv, cs%barotropic_CSp, cs%BT_cont%h_u, cs%BT_cont%h_v, &
539  obc=cs%OBC)
540  endif
541  if (showcalltree) call calltree_waypoint("done with continuity[BT_cont] (step_MOM_dyn_split_RK2)")
542  endif
543 
544  if (cs%BT_use_layer_fluxes) then
545  uh_ptr => uh_in ; vh_ptr => vh_in; u_ptr => u ; v_ptr => v
546  endif
547 
548  call cpu_clock_begin(id_clock_btstep)
549  if (calc_dtbt) call set_dtbt(g, gv, us, cs%barotropic_CSp, eta, cs%pbce)
550  if (showcalltree) call calltree_enter("btstep(), MOM_barotropic.F90")
551  ! This is the predictor step call to btstep.
552  call btstep(u, v, eta, dt, u_bc_accel, v_bc_accel, forces, cs%pbce, cs%eta_PF, &
553  u_av, v_av, cs%u_accel_bt, cs%v_accel_bt, eta_pred, cs%uhbt, cs%vhbt, &
554  g, gv, us, cs%barotropic_CSp, cs%visc_rem_u, cs%visc_rem_v, adp=cs%ADp, &
555  obc=cs%OBC, bt_cont=cs%BT_cont, eta_pf_start=eta_pf_start, &
556  taux_bot=taux_bot, tauy_bot=tauy_bot, &
557  uh0=uh_ptr, vh0=vh_ptr, u_uh0=u_ptr, v_vh0=v_ptr)
558  if (showcalltree) call calltree_leave("btstep()")
559  call cpu_clock_end(id_clock_btstep)
560 
561 ! up = u + dt_pred*( u_bc_accel + u_accel_bt )
562  dt_pred = dt * cs%be
563  call cpu_clock_begin(id_clock_mom_update)
564 
565  !$OMP parallel do default(shared)
566  do k=1,nz
567  do j=jsq,jeq ; do i=is,ie
568  vp(i,j,k) = g%mask2dCv(i,j) * (v(i,j,k) + dt_pred * &
569  (v_bc_accel(i,j,k) + cs%v_accel_bt(i,j,k)))
570  enddo ; enddo
571  do j=js,je ; do i=isq,ieq
572  up(i,j,k) = g%mask2dCu(i,j) * (u(i,j,k) + dt_pred * &
573  (u_bc_accel(i,j,k) + cs%u_accel_bt(i,j,k)))
574  enddo ; enddo
575  enddo
576  call cpu_clock_end(id_clock_mom_update)
577 
578  if (cs%debug) then
579  call uvchksum("Predictor 1 [uv]", up, vp, g%HI, haloshift=0, symmetric=sym, scale=us%L_T_to_m_s)
580  call hchksum(h, "Predictor 1 h", g%HI, haloshift=1, scale=gv%H_to_m)
581  call uvchksum("Predictor 1 [uv]h", uh, vh, g%HI,haloshift=2, &
582  symmetric=sym, scale=gv%H_to_m*us%L_to_m**2*us%s_to_T)
583 ! call MOM_state_chksum("Predictor 1", up, vp, h, uh, vh, G, GV, US, haloshift=1)
584  call mom_accel_chksum("Predictor accel", cs%CAu, cs%CAv, cs%PFu, cs%PFv, &
585  cs%diffu, cs%diffv, g, gv, us, cs%pbce, cs%u_accel_bt, cs%v_accel_bt, symmetric=sym)
586  call mom_state_chksum("Predictor 1 init", u, v, h, uh, vh, g, gv, us, haloshift=2, &
587  symmetric=sym)
588  call check_redundant("Predictor 1 up", up, vp, g)
589  call check_redundant("Predictor 1 uh", uh, vh, g)
590  endif
591 
592 ! up <- up + dt_pred d/dz visc d/dz up
593 ! u_av <- u_av + dt_pred d/dz visc d/dz u_av
594  call cpu_clock_begin(id_clock_vertvisc)
595  if (cs%debug) then
596  call uvchksum("0 before vertvisc: [uv]p", up, vp, g%HI,haloshift=0, symmetric=sym, scale=us%L_T_to_m_s)
597  endif
598  call vertvisc_coef(up, vp, h, forces, visc, dt_pred, g, gv, us, cs%vertvisc_CSp, &
599  cs%OBC)
600  call vertvisc(up, vp, h, forces, visc, dt_pred, cs%OBC, cs%ADp, cs%CDp, g, &
601  gv, us, cs%vertvisc_CSp, cs%taux_bot, cs%tauy_bot, waves=waves)
602  if (showcalltree) call calltree_waypoint("done with vertvisc (step_MOM_dyn_split_RK2)")
603  if (g%nonblocking_updates) then
604  call cpu_clock_end(id_clock_vertvisc)
605  call start_group_pass(cs%pass_uvp, g%Domain, clock=id_clock_pass)
606  call cpu_clock_begin(id_clock_vertvisc)
607  endif
608  call vertvisc_remnant(visc, cs%visc_rem_u, cs%visc_rem_v, dt_pred, g, gv, us, cs%vertvisc_CSp)
609  call cpu_clock_end(id_clock_vertvisc)
610 
611  call do_group_pass(cs%pass_visc_rem, g%Domain, clock=id_clock_pass)
612  if (g%nonblocking_updates) then
613  call complete_group_pass(cs%pass_uvp, g%Domain, clock=id_clock_pass)
614  else
615  call do_group_pass(cs%pass_uvp, g%Domain, clock=id_clock_pass)
616  endif
617 
618  ! uh = u_av * h
619  ! hp = h + dt * div . uh
620  call cpu_clock_begin(id_clock_continuity)
621  call continuity(up, vp, h, hp, uh, vh, dt, g, gv, us, cs%continuity_CSp, &
622  cs%uhbt, cs%vhbt, cs%OBC, cs%visc_rem_u, cs%visc_rem_v, &
623  u_av, v_av, bt_cont=cs%BT_cont)
624  call cpu_clock_end(id_clock_continuity)
625  if (showcalltree) call calltree_waypoint("done with continuity (step_MOM_dyn_split_RK2)")
626 
627  call do_group_pass(cs%pass_hp_uv, g%Domain, clock=id_clock_pass)
628 
629  if (associated(cs%OBC)) then
630 
631  if (cs%debug) &
632  call uvchksum("Pre OBC avg [uv]", u_av, v_av, g%HI, haloshift=1, symmetric=sym, scale=us%L_T_to_m_s)
633 
634  call radiation_open_bdry_conds(cs%OBC, u_av, u_old_rad_obc, v_av, v_old_rad_obc, g, us, dt_pred)
635 
636  if (cs%debug) &
637  call uvchksum("Post OBC avg [uv]", u_av, v_av, g%HI, haloshift=1, symmetric=sym, scale=us%L_T_to_m_s)
638 
639  ! These should be done with a pass that excludes uh & vh.
640 ! call do_group_pass(CS%pass_hp_uv, G%Domain, clock=id_clock_pass)
641  endif
642 
643  if (g%nonblocking_updates) then
644  call start_group_pass(cs%pass_av_uvh, g%Domain, clock=id_clock_pass)
645  endif
646 
647  ! h_av = (h + hp)/2
648  !$OMP parallel do default(shared)
649  do k=1,nz ; do j=js-2,je+2 ; do i=is-2,ie+2
650  h_av(i,j,k) = 0.5*(h(i,j,k) + hp(i,j,k))
651  enddo ; enddo ; enddo
652 
653  ! The correction phase of the time step starts here.
654  call enable_averages(dt, time_local, cs%diag)
655 
656  ! Calculate a revised estimate of the free-surface height correction to be
657  ! used in the next call to btstep. This call is at this point so that
658  ! hp can be changed if CS%begw /= 0.
659  ! eta_cor = ... (hidden inside CS%barotropic_CSp)
660  call cpu_clock_begin(id_clock_btcalc)
661  call bt_mass_source(hp, eta_pred, .false., g, gv, cs%barotropic_CSp)
662  call cpu_clock_end(id_clock_btcalc)
663 
664  if (cs%begw /= 0.0) then
665  ! hp <- (1-begw)*h_in + begw*hp
666  ! Back up hp to the value it would have had after a time-step of
667  ! begw*dt. hp is not used again until recalculated by continuity.
668  !$OMP parallel do default(shared)
669  do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1
670  hp(i,j,k) = (1.0-cs%begw)*h(i,j,k) + cs%begw*hp(i,j,k)
671  enddo ; enddo ; enddo
672 
673  ! PFu = d/dx M(hp,T,S)
674  ! pbce = dM/deta
675  call cpu_clock_begin(id_clock_pres)
676  call pressureforce(hp, tv, cs%PFu, cs%PFv, g, gv, us, cs%PressureForce_CSp, &
677  cs%ALE_CSp, p_surf, cs%pbce, cs%eta_PF)
678  call cpu_clock_end(id_clock_pres)
679  if (showcalltree) call calltree_waypoint("done with PressureForce[hp=(1-b).h+b.h] (step_MOM_dyn_split_RK2)")
680  endif
681 
682  if (g%nonblocking_updates) &
683  call complete_group_pass(cs%pass_av_uvh, g%Domain, clock=id_clock_pass)
684 
685  if (bt_cont_bt_thick) then
686  call btcalc(h, g, gv, cs%barotropic_CSp, cs%BT_cont%h_u, cs%BT_cont%h_v, &
687  obc=cs%OBC)
688  if (showcalltree) call calltree_waypoint("done with btcalc[BT_cont_BT_thick] (step_MOM_dyn_split_RK2)")
689  endif
690 
691  if (cs%debug) then
692  call mom_state_chksum("Predictor ", up, vp, hp, uh, vh, g, gv, us, symmetric=sym)
693  call uvchksum("Predictor avg [uv]", u_av, v_av, g%HI, haloshift=1, symmetric=sym, scale=us%L_T_to_m_s)
694  call hchksum(h_av, "Predictor avg h", g%HI, haloshift=0, scale=gv%H_to_m)
695  ! call MOM_state_chksum("Predictor avg ", u_av, v_av, h_av, uh, vh, G, GV, US)
696  call check_redundant("Predictor up ", up, vp, g)
697  call check_redundant("Predictor uh ", uh, vh, g)
698  endif
699 
700 ! diffu = horizontal viscosity terms (u_av)
701  call cpu_clock_begin(id_clock_horvisc)
702  call horizontal_viscosity(u_av, v_av, h_av, cs%diffu, cs%diffv, &
703  meke, varmix, g, gv, us, cs%hor_visc_CSp, &
704  obc=cs%OBC, bt=cs%barotropic_CSp, td=thickness_diffuse_csp, &
705  adp=cs%ADp)
706  call cpu_clock_end(id_clock_horvisc)
707  if (showcalltree) call calltree_waypoint("done with horizontal_viscosity (step_MOM_dyn_split_RK2)")
708 
709 ! CAu = -(f+zeta_av)/h_av vh + d/dx KE_av
710  call cpu_clock_begin(id_clock_cor)
711  call coradcalc(u_av, v_av, h_av, uh, vh, cs%CAu, cs%CAv, cs%OBC, cs%ADp, &
712  g, gv, us, cs%CoriolisAdv_CSp)
713  call cpu_clock_end(id_clock_cor)
714  if (showcalltree) call calltree_waypoint("done with CorAdCalc (step_MOM_dyn_split_RK2)")
715 
716 ! Calculate the momentum forcing terms for the barotropic equations.
717 
718 ! u_bc_accel = CAu + PFu + diffu(u[n-1])
719  call cpu_clock_begin(id_clock_btforce)
720  !$OMP parallel do default(shared)
721  do k=1,nz
722  do j=js,je ; do i=isq,ieq
723  u_bc_accel(i,j,k) = (cs%Cau(i,j,k) + cs%PFu(i,j,k)) + cs%diffu(i,j,k)
724  enddo ; enddo
725  do j=jsq,jeq ; do i=is,ie
726  v_bc_accel(i,j,k) = (cs%Cav(i,j,k) + cs%PFv(i,j,k)) + cs%diffv(i,j,k)
727  enddo ; enddo
728  enddo
729  if (associated(cs%OBC)) then
730  call open_boundary_zero_normal_flow(cs%OBC, g, u_bc_accel, v_bc_accel)
731  endif
732  call cpu_clock_end(id_clock_btforce)
733 
734  if (cs%debug) then
735  call mom_accel_chksum("corr pre-btstep accel", cs%CAu, cs%CAv, cs%PFu, cs%PFv, &
736  cs%diffu, cs%diffv, g, gv, us, cs%pbce, u_bc_accel, v_bc_accel, &
737  symmetric=sym)
738  call check_redundant("corr pre-btstep CS%Ca ", cs%Cau, cs%Cav, g)
739  call check_redundant("corr pre-btstep CS%PF ", cs%PFu, cs%PFv, g)
740  call check_redundant("corr pre-btstep CS%diff ", cs%diffu, cs%diffv, g)
741  call check_redundant("corr pre-btstep u_bc_accel ", u_bc_accel, v_bc_accel, g)
742  endif
743 
744  ! u_accel_bt = layer accelerations due to barotropic solver
745  ! pbce = dM/deta
746  call cpu_clock_begin(id_clock_btstep)
747  if (cs%BT_use_layer_fluxes) then
748  uh_ptr => uh ; vh_ptr => vh ; u_ptr => u_av ; v_ptr => v_av
749  endif
750 
751  if (showcalltree) call calltree_enter("btstep(), MOM_barotropic.F90")
752  ! This is the corrector step call to btstep.
753  call btstep(u, v, eta, dt, u_bc_accel, v_bc_accel, forces, cs%pbce, &
754  cs%eta_PF, u_av, v_av, cs%u_accel_bt, cs%v_accel_bt, &
755  eta_pred, cs%uhbt, cs%vhbt, g, gv, us, cs%barotropic_CSp, &
756  cs%visc_rem_u, cs%visc_rem_v, etaav=eta_av, adp=cs%ADp, &
757  obc=cs%OBC, bt_cont = cs%BT_cont, eta_pf_start=eta_pf_start, &
758  taux_bot=taux_bot, tauy_bot=tauy_bot, &
759  uh0=uh_ptr, vh0=vh_ptr, u_uh0=u_ptr, v_vh0=v_ptr)
760  do j=js,je ; do i=is,ie ; eta(i,j) = eta_pred(i,j) ; enddo ; enddo
761 
762  call cpu_clock_end(id_clock_btstep)
763  if (showcalltree) call calltree_leave("btstep()")
764 
765  if (cs%debug) then
766  call check_redundant("u_accel_bt ", cs%u_accel_bt, cs%v_accel_bt, g)
767  endif
768 
769  ! u = u + dt*( u_bc_accel + u_accel_bt )
770  call cpu_clock_begin(id_clock_mom_update)
771  !$OMP parallel do default(shared)
772  do k=1,nz
773  do j=js,je ; do i=isq,ieq
774  u(i,j,k) = g%mask2dCu(i,j) * (u(i,j,k) + dt * &
775  (u_bc_accel(i,j,k) + cs%u_accel_bt(i,j,k)))
776  enddo ; enddo
777  do j=jsq,jeq ; do i=is,ie
778  v(i,j,k) = g%mask2dCv(i,j) * (v(i,j,k) + dt * &
779  (v_bc_accel(i,j,k) + cs%v_accel_bt(i,j,k)))
780  enddo ; enddo
781  enddo
782  call cpu_clock_end(id_clock_mom_update)
783 
784  if (cs%debug) then
785  call uvchksum("Corrector 1 [uv]", u, v, g%HI,haloshift=0, symmetric=sym, scale=us%L_T_to_m_s)
786  call hchksum(h, "Corrector 1 h", g%HI, haloshift=2, scale=gv%H_to_m)
787  call uvchksum("Corrector 1 [uv]h", uh, vh, g%HI, haloshift=2, &
788  symmetric=sym, scale=gv%H_to_m*us%L_to_m**2*us%s_to_T)
789  ! call MOM_state_chksum("Corrector 1", u, v, h, uh, vh, G, GV, US, haloshift=1)
790  call mom_accel_chksum("Corrector accel", cs%CAu, cs%CAv, cs%PFu, cs%PFv, &
791  cs%diffu, cs%diffv, g, gv, us, cs%pbce, cs%u_accel_bt, cs%v_accel_bt, &
792  symmetric=sym)
793  endif
794 
795  ! u <- u + dt d/dz visc d/dz u
796  ! u_av <- u_av + dt d/dz visc d/dz u_av
797  call cpu_clock_begin(id_clock_vertvisc)
798  call vertvisc_coef(u, v, h, forces, visc, dt, g, gv, us, cs%vertvisc_CSp, cs%OBC)
799  call vertvisc(u, v, h, forces, visc, dt, cs%OBC, cs%ADp, cs%CDp, g, gv, us, &
800  cs%vertvisc_CSp, cs%taux_bot, cs%tauy_bot,waves=waves)
801  if (g%nonblocking_updates) then
802  call cpu_clock_end(id_clock_vertvisc)
803  call start_group_pass(cs%pass_uv, g%Domain, clock=id_clock_pass)
804  call cpu_clock_begin(id_clock_vertvisc)
805  endif
806  call vertvisc_remnant(visc, cs%visc_rem_u, cs%visc_rem_v, dt, g, gv, us, cs%vertvisc_CSp)
807  call cpu_clock_end(id_clock_vertvisc)
808  if (showcalltree) call calltree_waypoint("done with vertvisc (step_MOM_dyn_split_RK2)")
809 
810 ! Later, h_av = (h_in + h_out)/2, but for now use h_av to store h_in.
811  !$OMP parallel do default(shared)
812  do k=1,nz ; do j=js-2,je+2 ; do i=is-2,ie+2
813  h_av(i,j,k) = h(i,j,k)
814  enddo ; enddo ; enddo
815 
816  call do_group_pass(cs%pass_visc_rem, g%Domain, clock=id_clock_pass)
817  if (g%nonblocking_updates) then
818  call complete_group_pass(cs%pass_uv, g%Domain, clock=id_clock_pass)
819  else
820  call do_group_pass(cs%pass_uv, g%Domain, clock=id_clock_pass)
821  endif
822 
823  ! uh = u_av * h
824  ! h = h + dt * div . uh
825  ! u_av and v_av adjusted so their mass transports match uhbt and vhbt.
826  call cpu_clock_begin(id_clock_continuity)
827  call continuity(u, v, h, h, uh, vh, dt, g, gv, us, cs%continuity_CSp, &
828  cs%uhbt, cs%vhbt, cs%OBC, cs%visc_rem_u, cs%visc_rem_v, u_av, v_av)
829  call cpu_clock_end(id_clock_continuity)
830  call do_group_pass(cs%pass_h, g%Domain, clock=id_clock_pass)
831  ! Whenever thickness changes let the diag manager know, target grids
832  ! for vertical remapping may need to be regenerated.
833  call diag_update_remap_grids(cs%diag)
834  if (showcalltree) call calltree_waypoint("done with continuity (step_MOM_dyn_split_RK2)")
835 
836  if (g%nonblocking_updates) then
837  call start_group_pass(cs%pass_av_uvh, g%Domain, clock=id_clock_pass)
838  else
839  call do_group_pass(cs%pass_av_uvh, g%domain, clock=id_clock_pass)
840  endif
841 
842  if (associated(cs%OBC)) then
843  call radiation_open_bdry_conds(cs%OBC, u, u_old_rad_obc, v, v_old_rad_obc, g, us, dt)
844  endif
845 
846 ! h_av = (h_in + h_out)/2 . Going in to this line, h_av = h_in.
847  !$OMP parallel do default(shared)
848  do k=1,nz ; do j=js-2,je+2 ; do i=is-2,ie+2
849  h_av(i,j,k) = 0.5*(h_av(i,j,k) + h(i,j,k))
850  enddo ; enddo ; enddo
851 
852  if (g%nonblocking_updates) &
853  call complete_group_pass(cs%pass_av_uvh, g%Domain, clock=id_clock_pass)
854 
855  !$OMP parallel do default(shared)
856  do k=1,nz
857  do j=js-2,je+2 ; do i=isq-2,ieq+2
858  uhtr(i,j,k) = uhtr(i,j,k) + uh(i,j,k)*dt
859  enddo ; enddo
860  do j=jsq-2,jeq+2 ; do i=is-2,ie+2
861  vhtr(i,j,k) = vhtr(i,j,k) + vh(i,j,k)*dt
862  enddo ; enddo
863  enddo
864 
865  ! The time-averaged free surface height has already been set by the last
866  ! call to btstep.
867 
868  ! Here various terms used in to update the momentum equations are
869  ! offered for time averaging.
870  if (cs%id_PFu > 0) call post_data(cs%id_PFu, cs%PFu, cs%diag)
871  if (cs%id_PFv > 0) call post_data(cs%id_PFv, cs%PFv, cs%diag)
872  if (cs%id_CAu > 0) call post_data(cs%id_CAu, cs%CAu, cs%diag)
873  if (cs%id_CAv > 0) call post_data(cs%id_CAv, cs%CAv, cs%diag)
874 
875  ! Here the thickness fluxes are offered for time averaging.
876  if (cs%id_uh > 0) call post_data(cs%id_uh , uh, cs%diag)
877  if (cs%id_vh > 0) call post_data(cs%id_vh , vh, cs%diag)
878  if (cs%id_uav > 0) call post_data(cs%id_uav, u_av, cs%diag)
879  if (cs%id_vav > 0) call post_data(cs%id_vav, v_av, cs%diag)
880  if (cs%id_u_BT_accel > 0) call post_data(cs%id_u_BT_accel, cs%u_accel_bt, cs%diag)
881  if (cs%id_v_BT_accel > 0) call post_data(cs%id_v_BT_accel, cs%v_accel_bt, cs%diag)
882 
883  ! Diagnostics for terms multiplied by fractional thicknesses
884 
885  ! 3D diagnostics hf_PFu etc. are commented because there is no clarity on proper remapping grid option.
886  ! The code is retained for degugging purposes in the future.
887  !if (CS%id_hf_PFu > 0) then
888  ! allocate(hf_PFu(G%IsdB:G%IedB,G%jsd:G%jed,G%ke))
889  ! do k=1,nz ; do j=js,je ; do I=Isq,Ieq
890  ! hf_PFu(I,j,k) = CS%PFu(I,j,k) * CS%ADp%diag_hfrac_u(I,j,k)
891  ! enddo ; enddo ; enddo
892  ! call post_data(CS%id_hf_PFu, hf_PFu, CS%diag)
893  !endif
894  !if (CS%id_hf_PFv > 0) then
895  ! allocate(hf_PFv(G%isd:G%ied,G%JsdB:G%JedB,G%ke))
896  ! do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
897  ! hf_PFv(i,J,k) = CS%PFv(i,J,k) * CS%ADp%diag_hfrac_v(i,J,k)
898  ! enddo ; enddo ; enddo
899  ! call post_data(CS%id_hf_PFv, hf_PFv, CS%diag)
900  !endif
901  if (cs%id_hf_PFu_2d > 0) then
902  allocate(hf_pfu_2d(g%IsdB:g%IedB,g%jsd:g%jed))
903  hf_pfu_2d(:,:) = 0.0
904  do k=1,nz ; do j=js,je ; do i=isq,ieq
905  hf_pfu_2d(i,j) = hf_pfu_2d(i,j) + cs%PFu(i,j,k) * cs%ADp%diag_hfrac_u(i,j,k)
906  enddo ; enddo ; enddo
907  call post_data(cs%id_hf_PFu_2d, hf_pfu_2d, cs%diag)
908  deallocate(hf_pfu_2d)
909  endif
910  if (cs%id_hf_PFv_2d > 0) then
911  allocate(hf_pfv_2d(g%isd:g%ied,g%JsdB:g%JedB))
912  hf_pfv_2d(:,:) = 0.0
913  do k=1,nz ; do j=jsq,jeq ; do i=is,ie
914  hf_pfv_2d(i,j) = hf_pfv_2d(i,j) + cs%PFv(i,j,k) * cs%ADp%diag_hfrac_v(i,j,k)
915  enddo ; enddo ; enddo
916  call post_data(cs%id_hf_PFv_2d, hf_pfv_2d, cs%diag)
917  deallocate(hf_pfv_2d)
918  endif
919 
920  !if (CS%id_hf_CAu > 0) then
921  ! allocate(hf_CAu(G%IsdB:G%IedB,G%jsd:G%jed,G%ke))
922  ! do k=1,nz ; do j=js,je ; do I=Isq,Ieq
923  ! hf_CAu(I,j,k) = CS%CAu(I,j,k) * CS%ADp%diag_hfrac_u(I,j,k)
924  ! enddo ; enddo ; enddo
925  ! call post_data(CS%id_hf_CAu, hf_CAu, CS%diag)
926  !endif
927  !if (CS%id_hf_CAv > 0) then
928  ! allocate(hf_CAv(G%isd:G%ied,G%JsdB:G%JedB,G%ke))
929  ! do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
930  ! hf_CAv(i,J,k) = CS%CAv(i,J,k) * CS%ADp%diag_hfrac_v(i,J,k)
931  ! enddo ; enddo ; enddo
932  ! call post_data(CS%id_hf_CAv, hf_CAv, CS%diag)
933  !endif
934  if (cs%id_hf_CAu_2d > 0) then
935  allocate(hf_cau_2d(g%IsdB:g%IedB,g%jsd:g%jed))
936  hf_cau_2d(:,:) = 0.0
937  do k=1,nz ; do j=js,je ; do i=isq,ieq
938  hf_cau_2d(i,j) = hf_cau_2d(i,j) + cs%CAu(i,j,k) * cs%ADp%diag_hfrac_u(i,j,k)
939  enddo ; enddo ; enddo
940  call post_data(cs%id_hf_CAu_2d, hf_cau_2d, cs%diag)
941  deallocate(hf_cau_2d)
942  endif
943  if (cs%id_hf_CAv_2d > 0) then
944  allocate(hf_cav_2d(g%isd:g%ied,g%JsdB:g%JedB))
945  hf_cav_2d(:,:) = 0.0
946  do k=1,nz ; do j=jsq,jeq ; do i=is,ie
947  hf_cav_2d(i,j) = hf_cav_2d(i,j) + cs%CAv(i,j,k) * cs%ADp%diag_hfrac_v(i,j,k)
948  enddo ; enddo ; enddo
949  call post_data(cs%id_hf_CAv_2d, hf_cav_2d, cs%diag)
950  deallocate(hf_cav_2d)
951  endif
952 
953  !if (CS%id_hf_u_BT_accel > 0) then
954  ! allocate(hf_u_BT_accel(G%IsdB:G%IedB,G%jsd:G%jed,G%ke))
955  ! do k=1,nz ; do j=js,je ; do I=Isq,Ieq
956  ! hf_u_BT_accel(I,j,k) = CS%u_accel_bt(I,j,k) * CS%ADp%diag_hfrac_u(I,j,k)
957  ! enddo ; enddo ; enddo
958  ! call post_data(CS%id_hf_u_BT_accel, hf_u_BT_accel, CS%diag)
959  !endif
960  !if (CS%id_hf_v_BT_accel > 0) then
961  ! allocate(hf_v_BT_accel(G%isd:G%ied,G%JsdB:G%JedB,G%ke))
962  ! do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
963  ! hf_v_BT_accel(i,J,k) = CS%v_accel_bt(i,J,k) * CS%ADp%diag_hfrac_v(i,J,k)
964  ! enddo ; enddo ; enddo
965  ! call post_data(CS%id_hf_v_BT_accel, hf_v_BT_accel, CS%diag)
966  !endif
967  if (cs%id_hf_u_BT_accel_2d > 0) then
968  allocate(hf_u_bt_accel_2d(g%IsdB:g%IedB,g%jsd:g%jed))
969  hf_u_bt_accel_2d(:,:) = 0.0
970  do k=1,nz ; do j=js,je ; do i=isq,ieq
971  hf_u_bt_accel_2d(i,j) = hf_u_bt_accel_2d(i,j) + cs%u_accel_bt(i,j,k) * cs%ADp%diag_hfrac_u(i,j,k)
972  enddo ; enddo ; enddo
973  call post_data(cs%id_hf_u_BT_accel_2d, hf_u_bt_accel_2d, cs%diag)
974  deallocate(hf_u_bt_accel_2d)
975  endif
976  if (cs%id_hf_v_BT_accel_2d > 0) then
977  allocate(hf_v_bt_accel_2d(g%isd:g%ied,g%JsdB:g%JedB))
978  hf_v_bt_accel_2d(:,:) = 0.0
979  do k=1,nz ; do j=jsq,jeq ; do i=is,ie
980  hf_v_bt_accel_2d(i,j) = hf_v_bt_accel_2d(i,j) + cs%v_accel_bt(i,j,k) * cs%ADp%diag_hfrac_v(i,j,k)
981  enddo ; enddo ; enddo
982  call post_data(cs%id_hf_v_BT_accel_2d, hf_v_bt_accel_2d, cs%diag)
983  deallocate(hf_v_bt_accel_2d)
984  endif
985 
986  if (cs%debug) then
987  call mom_state_chksum("Corrector ", u, v, h, uh, vh, g, gv, us, symmetric=sym)
988  call uvchksum("Corrector avg [uv]", u_av, v_av, g%HI, haloshift=1, symmetric=sym, scale=us%L_T_to_m_s)
989  call hchksum(h_av, "Corrector avg h", g%HI, haloshift=1, scale=gv%H_to_m)
990  ! call MOM_state_chksum("Corrector avg ", u_av, v_av, h_av, uh, vh, G, GV, US)
991  endif
992 
993  if (showcalltree) call calltree_leave("step_MOM_dyn_split_RK2()")
994 
995 end subroutine step_mom_dyn_split_rk2
996 
997 !> This subroutine sets up any auxiliary restart variables that are specific
998 !! to the unsplit time stepping scheme. All variables registered here should
999 !! have the ability to be recreated if they are not present in a restart file.
1000 subroutine register_restarts_dyn_split_rk2(HI, GV, param_file, CS, restart_CS, uh, vh)
1001  type(hor_index_type), intent(in) :: hi !< Horizontal index structure
1002  type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
1003  type(param_file_type), intent(in) :: param_file !< parameter file
1004  type(mom_dyn_split_rk2_cs), pointer :: cs !< module control structure
1005  type(mom_restart_cs), pointer :: restart_cs !< restart control structure
1006  real, dimension(SZIB_(HI),SZJ_(HI),SZK_(GV)), &
1007  target, intent(inout) :: uh !< zonal volume/mass transport [H L2 T-1 ~> m3 s-1 or kg s-1]
1008  real, dimension(SZI_(HI),SZJB_(HI),SZK_(GV)), &
1009  target, intent(inout) :: vh !< merid volume/mass transport [H L2 T-1 ~> m3 s-1 or kg s-1]
1010 
1011  type(vardesc) :: vd(2)
1012  character(len=40) :: mdl = "MOM_dynamics_split_RK2" ! This module's name.
1013  character(len=48) :: thickness_units, flux_units
1014 
1015  integer :: isd, ied, jsd, jed, nz, isdb, iedb, jsdb, jedb
1016 
1017  isd = hi%isd ; ied = hi%ied ; jsd = hi%jsd ; jed = hi%jed ; nz = gv%ke
1018  isdb = hi%IsdB ; iedb = hi%IedB ; jsdb = hi%JsdB ; jedb = hi%JedB
1019 
1020  ! This is where a control structure specific to this module would be allocated.
1021  if (associated(cs)) then
1022  call mom_error(warning, "register_restarts_dyn_split_RK2 called with an associated "// &
1023  "control structure.")
1024  return
1025  endif
1026  allocate(cs)
1027 
1028  alloc_(cs%diffu(isdb:iedb,jsd:jed,nz)) ; cs%diffu(:,:,:) = 0.0
1029  alloc_(cs%diffv(isd:ied,jsdb:jedb,nz)) ; cs%diffv(:,:,:) = 0.0
1030  alloc_(cs%CAu(isdb:iedb,jsd:jed,nz)) ; cs%CAu(:,:,:) = 0.0
1031  alloc_(cs%CAv(isd:ied,jsdb:jedb,nz)) ; cs%CAv(:,:,:) = 0.0
1032  alloc_(cs%PFu(isdb:iedb,jsd:jed,nz)) ; cs%PFu(:,:,:) = 0.0
1033  alloc_(cs%PFv(isd:ied,jsdb:jedb,nz)) ; cs%PFv(:,:,:) = 0.0
1034 
1035  alloc_(cs%eta(isd:ied,jsd:jed)) ; cs%eta(:,:) = 0.0
1036  alloc_(cs%u_av(isdb:iedb,jsd:jed,nz)) ; cs%u_av(:,:,:) = 0.0
1037  alloc_(cs%v_av(isd:ied,jsdb:jedb,nz)) ; cs%v_av(:,:,:) = 0.0
1038  alloc_(cs%h_av(isd:ied,jsd:jed,nz)) ; cs%h_av(:,:,:) = gv%Angstrom_H
1039 
1040  thickness_units = get_thickness_units(gv)
1041  flux_units = get_flux_units(gv)
1042 
1043  if (gv%Boussinesq) then
1044  vd(1) = var_desc("sfc",thickness_units,"Free surface Height",'h','1')
1045  else
1046  vd(1) = var_desc("p_bot",thickness_units,"Bottom Pressure",'h','1')
1047  endif
1048  call register_restart_field(cs%eta, vd(1), .false., restart_cs)
1049 
1050  vd(1) = var_desc("u2","m s-1","Auxiliary Zonal velocity",'u','L')
1051  vd(2) = var_desc("v2","m s-1","Auxiliary Meridional velocity",'v','L')
1052  call register_restart_pair(cs%u_av, cs%v_av, vd(1), vd(2), .false., restart_cs)
1053 
1054  vd(1) = var_desc("h2",thickness_units,"Auxiliary Layer Thickness",'h','L')
1055  call register_restart_field(cs%h_av, vd(1), .false., restart_cs)
1056 
1057  vd(1) = var_desc("uh",flux_units,"Zonal thickness flux",'u','L')
1058  vd(2) = var_desc("vh",flux_units,"Meridional thickness flux",'v','L')
1059  call register_restart_pair(uh, vh, vd(1), vd(2), .false., restart_cs)
1060 
1061  vd(1) = var_desc("diffu","m s-2","Zonal horizontal viscous acceleration",'u','L')
1062  vd(2) = var_desc("diffv","m s-2","Meridional horizontal viscous acceleration",'v','L')
1063  call register_restart_pair(cs%diffu, cs%diffv, vd(1), vd(2), .false., restart_cs)
1064 
1065  call register_barotropic_restarts(hi, gv, param_file, cs%barotropic_CSp, &
1066  restart_cs)
1067 
1068 end subroutine register_restarts_dyn_split_rk2
1069 
1070 !> This subroutine initializes all of the variables that are used by this
1071 !! dynamic core, including diagnostics and the cpu clocks.
1072 subroutine initialize_dyn_split_rk2(u, v, h, uh, vh, eta, Time, G, GV, US, param_file, &
1073  diag, CS, restart_CS, dt, Accel_diag, Cont_diag, MIS, &
1074  VarMix, MEKE, thickness_diffuse_CSp, &
1075  OBC, update_OBC_CSp, ALE_CSp, setVisc_CSp, &
1076  visc, dirs, ntrunc, calc_dtbt, cont_stencil)
1077  type(ocean_grid_type), intent(inout) :: g !< ocean grid structure
1078  type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
1079  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1080  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
1081  intent(inout) :: u !< zonal velocity [L T-1 ~> m s-1]
1082  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1083  intent(inout) :: v !< merid velocity [L T-1 ~> m s-1]
1084  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) , intent(inout) :: h !< layer thickness [H ~> m or kg m-2]
1085  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
1086  target, intent(inout) :: uh !< zonal volume/mass transport [H L2 T-1 ~> m3 s-1 or kg s-1]
1087  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1088  target, intent(inout) :: vh !< merid volume/mass transport [H L2 T-1 ~> m3 s-1 or kg s-1]
1089  real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: eta !< free surface height or column mass [H ~> m or kg m-2]
1090  type(time_type), target, intent(in) :: time !< current model time
1091  type(param_file_type), intent(in) :: param_file !< parameter file for parsing
1092  type(diag_ctrl), target, intent(inout) :: diag !< to control diagnostics
1093  type(mom_dyn_split_rk2_cs), pointer :: cs !< module control structure
1094  type(mom_restart_cs), pointer :: restart_cs !< restart control structure
1095  real, intent(in) :: dt !< time step [T ~> s]
1096  type(accel_diag_ptrs), target, intent(inout) :: accel_diag !< points to momentum equation terms for
1097  !! budget analysis
1098  type(cont_diag_ptrs), target, intent(inout) :: cont_diag !< points to terms in continuity equation
1099  type(ocean_internal_state), intent(inout) :: mis !< "MOM6 internal state" used to pass
1100  !! diagnostic pointers
1101  type(varmix_cs), pointer :: varmix !< points to spatially variable viscosities
1102  type(meke_type), pointer :: meke !< points to mesoscale eddy kinetic energy fields
1103 ! type(Barotropic_CS), pointer :: Barotropic_CSp !< Pointer to the control structure for
1104 ! !! the barotropic module
1105  type(thickness_diffuse_cs), pointer :: thickness_diffuse_csp !< Pointer to the control structure
1106  !! used for the isopycnal height diffusive transport.
1107  type(ocean_obc_type), pointer :: obc !< points to OBC related fields
1108  type(update_obc_cs), pointer :: update_obc_csp !< points to OBC update related fields
1109  type(ale_cs), pointer :: ale_csp !< points to ALE control structure
1110  type(set_visc_cs), pointer :: setvisc_csp !< points to the set_visc control structure.
1111  type(vertvisc_type), intent(inout) :: visc !< vertical viscosities, bottom drag, and related
1112  type(directories), intent(in) :: dirs !< contains directory paths
1113  integer, target, intent(inout) :: ntrunc !< A target for the variable that records
1114  !! the number of times the velocity is
1115  !! truncated (this should be 0).
1116  logical, intent(out) :: calc_dtbt !< If true, recalculate the barotropic time step
1117  integer, optional, intent(out) :: cont_stencil !< The stencil for thickness
1118  !! from the continuity solver.
1119 
1120  ! local variables
1121  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_tmp
1122  character(len=40) :: mdl = "MOM_dynamics_split_RK2" ! This module's name.
1123  ! This include declares and sets the variable "version".
1124 # include "version_variable.h"
1125  character(len=48) :: thickness_units, flux_units, eta_rest_name
1126  real :: h_rescale ! A rescaling factor for thicknesses from the representation in
1127  ! a restart file to the internal representation in this run.
1128  real :: vel_rescale ! A rescaling factor for velocities from the representation in
1129  ! a restart file to the internal representation in this run.
1130  real :: uh_rescale ! A rescaling factor for thickness transports from the representation in
1131  ! a restart file to the internal representation in this run.
1132  real :: accel_rescale ! A rescaling factor for accelerations from the representation in
1133  ! a restart file to the internal representation in this run.
1134  real :: h_convert
1135  type(group_pass_type) :: pass_av_h_uvh
1136  logical :: use_tides, debug_truncations
1137 
1138  integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz
1139  integer :: isdb, iedb, jsdb, jedb
1140  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
1141  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
1142  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
1143 
1144  if (.not.associated(cs)) call mom_error(fatal, &
1145  "initialize_dyn_split_RK2 called with an unassociated control structure.")
1146  if (cs%module_is_initialized) then
1147  call mom_error(warning, "initialize_dyn_split_RK2 called with a control "// &
1148  "structure that has already been initialized.")
1149  return
1150  endif
1151  cs%module_is_initialized = .true.
1152 
1153  cs%diag => diag
1154 
1155  call log_version(param_file, mdl, version, "")
1156  call get_param(param_file, mdl, "TIDES", use_tides, &
1157  "If true, apply tidal momentum forcing.", default=.false.)
1158  call get_param(param_file, mdl, "BE", cs%be, &
1159  "If SPLIT is true, BE determines the relative weighting "//&
1160  "of a 2nd-order Runga-Kutta baroclinic time stepping "//&
1161  "scheme (0.5) and a backward Euler scheme (1) that is "//&
1162  "used for the Coriolis and inertial terms. BE may be "//&
1163  "from 0.5 to 1, but instability may occur near 0.5. "//&
1164  "BE is also applicable if SPLIT is false and USE_RK2 "//&
1165  "is true.", units="nondim", default=0.6)
1166  call get_param(param_file, mdl, "BEGW", cs%begw, &
1167  "If SPLIT is true, BEGW is a number from 0 to 1 that "//&
1168  "controls the extent to which the treatment of gravity "//&
1169  "waves is forward-backward (0) or simulated backward "//&
1170  "Euler (1). 0 is almost always used. "//&
1171  "If SPLIT is false and USE_RK2 is true, BEGW can be "//&
1172  "between 0 and 0.5 to damp gravity waves.", &
1173  units="nondim", default=0.0)
1174 
1175  call get_param(param_file, mdl, "SPLIT_BOTTOM_STRESS", cs%split_bottom_stress, &
1176  "If true, provide the bottom stress calculated by the "//&
1177  "vertical viscosity to the barotropic solver.", default=.false.)
1178  call get_param(param_file, mdl, "BT_USE_LAYER_FLUXES", cs%BT_use_layer_fluxes, &
1179  "If true, use the summed layered fluxes plus an "//&
1180  "adjustment due to the change in the barotropic velocity "//&
1181  "in the barotropic continuity equation.", default=.true.)
1182  call get_param(param_file, mdl, "DEBUG", cs%debug, &
1183  "If true, write out verbose debugging data.", &
1184  default=.false., debuggingparam=.true.)
1185  call get_param(param_file, mdl, "DEBUG_OBC", cs%debug_OBC, default=.false.)
1186  call get_param(param_file, mdl, "DEBUG_TRUNCATIONS", debug_truncations, &
1187  default=.false.)
1188 
1189  allocate(cs%taux_bot(isdb:iedb,jsd:jed)) ; cs%taux_bot(:,:) = 0.0
1190  allocate(cs%tauy_bot(isd:ied,jsdb:jedb)) ; cs%tauy_bot(:,:) = 0.0
1191 
1192  alloc_(cs%uhbt(isdb:iedb,jsd:jed)) ; cs%uhbt(:,:) = 0.0
1193  alloc_(cs%vhbt(isd:ied,jsdb:jedb)) ; cs%vhbt(:,:) = 0.0
1194  alloc_(cs%visc_rem_u(isdb:iedb,jsd:jed,nz)) ; cs%visc_rem_u(:,:,:) = 0.0
1195  alloc_(cs%visc_rem_v(isd:ied,jsdb:jedb,nz)) ; cs%visc_rem_v(:,:,:) = 0.0
1196  alloc_(cs%eta_PF(isd:ied,jsd:jed)) ; cs%eta_PF(:,:) = 0.0
1197  alloc_(cs%pbce(isd:ied,jsd:jed,nz)) ; cs%pbce(:,:,:) = 0.0
1198 
1199  alloc_(cs%u_accel_bt(isdb:iedb,jsd:jed,nz)) ; cs%u_accel_bt(:,:,:) = 0.0
1200  alloc_(cs%v_accel_bt(isd:ied,jsdb:jedb,nz)) ; cs%v_accel_bt(:,:,:) = 0.0
1201 
1202  mis%diffu => cs%diffu
1203  mis%diffv => cs%diffv
1204  mis%PFu => cs%PFu
1205  mis%PFv => cs%PFv
1206  mis%CAu => cs%CAu
1207  mis%CAv => cs%CAv
1208  mis%pbce => cs%pbce
1209  mis%u_accel_bt => cs%u_accel_bt
1210  mis%v_accel_bt => cs%v_accel_bt
1211  mis%u_av => cs%u_av
1212  mis%v_av => cs%v_av
1213 
1214  cs%ADp => accel_diag
1215  cs%CDp => cont_diag
1216  accel_diag%diffu => cs%diffu
1217  accel_diag%diffv => cs%diffv
1218  accel_diag%PFu => cs%PFu
1219  accel_diag%PFv => cs%PFv
1220  accel_diag%CAu => cs%CAu
1221  accel_diag%CAv => cs%CAv
1222 
1223 ! Accel_diag%pbce => CS%pbce
1224 ! Accel_diag%u_accel_bt => CS%u_accel_bt ; Accel_diag%v_accel_bt => CS%v_accel_bt
1225 ! Accel_diag%u_av => CS%u_av ; Accel_diag%v_av => CS%v_av
1226 
1227  id_clock_pass_init = cpu_clock_id('(Ocean init message passing)', &
1228  grain=clock_routine)
1229 
1230  call continuity_init(time, g, gv, us, param_file, diag, cs%continuity_CSp)
1231  if (present(cont_stencil)) cont_stencil = continuity_stencil(cs%continuity_CSp)
1232  call coriolisadv_init(time, g, gv, us, param_file, diag, cs%ADp, cs%CoriolisAdv_CSp)
1233  if (use_tides) call tidal_forcing_init(time, g, param_file, cs%tides_CSp)
1234  call pressureforce_init(time, g, gv, us, param_file, diag, cs%PressureForce_CSp, &
1235  cs%tides_CSp)
1236  call hor_visc_init(time, g, us, param_file, diag, cs%hor_visc_CSp, meke, adp=cs%ADp)
1237  call vertvisc_init(mis, time, g, gv, us, param_file, diag, cs%ADp, dirs, &
1238  ntrunc, cs%vertvisc_CSp)
1239  if (.not.associated(setvisc_csp)) call mom_error(fatal, &
1240  "initialize_dyn_split_RK2 called with setVisc_CSp unassociated.")
1241  cs%set_visc_CSp => setvisc_csp
1242  call updatecfltruncationvalue(time, cs%vertvisc_CSp, &
1243  activate=is_new_run(restart_cs) )
1244 
1245  if (associated(ale_csp)) cs%ALE_CSp => ale_csp
1246  if (associated(obc)) then
1247  cs%OBC => obc
1248  if (obc%ramp) call update_obc_ramp(time, cs%OBC, &
1249  activate=is_new_run(restart_cs) )
1250  endif
1251  if (associated(update_obc_csp)) cs%update_OBC_CSp => update_obc_csp
1252 
1253  eta_rest_name = "sfc" ; if (.not.gv%Boussinesq) eta_rest_name = "p_bot"
1254  if (.not. query_initialized(cs%eta, trim(eta_rest_name), restart_cs)) then
1255  ! Estimate eta based on the layer thicknesses - h. With the Boussinesq
1256  ! approximation, eta is the free surface height anomaly, while without it
1257  ! eta is the mass of ocean per unit area. eta always has the same
1258  ! dimensions as h, either m or kg m-3.
1259  ! CS%eta(:,:) = 0.0 already from initialization.
1260  if (gv%Boussinesq) then
1261  do j=js,je ; do i=is,ie ; cs%eta(i,j) = -gv%Z_to_H * g%bathyT(i,j) ; enddo ; enddo
1262  endif
1263  do k=1,nz ; do j=js,je ; do i=is,ie
1264  cs%eta(i,j) = cs%eta(i,j) + h(i,j,k)
1265  enddo ; enddo ; enddo
1266  elseif ((gv%m_to_H_restart /= 0.0) .and. (gv%m_to_H_restart /= gv%m_to_H)) then
1267  h_rescale = gv%m_to_H / gv%m_to_H_restart
1268  do j=js,je ; do i=is,ie ; cs%eta(i,j) = h_rescale * cs%eta(i,j) ; enddo ; enddo
1269  endif
1270  ! Copy eta into an output array.
1271  do j=js,je ; do i=is,ie ; eta(i,j) = cs%eta(i,j) ; enddo ; enddo
1272 
1273  call barotropic_init(u, v, h, cs%eta, time, g, gv, us, param_file, diag, &
1274  cs%barotropic_CSp, restart_cs, calc_dtbt, cs%BT_cont, &
1275  cs%tides_CSp)
1276 
1277  if (.not. query_initialized(cs%diffu,"diffu",restart_cs) .or. &
1278  .not. query_initialized(cs%diffv,"diffv",restart_cs)) then
1279  call horizontal_viscosity(u, v, h, cs%diffu, cs%diffv, meke, varmix, &
1280  g, gv, us, cs%hor_visc_CSp, &
1281  obc=cs%OBC, bt=cs%barotropic_CSp, &
1282  td=thickness_diffuse_csp)
1283  else
1284  if ( (us%s_to_T_restart * us%m_to_L_restart /= 0.0) .and. &
1285  (us%m_to_L * us%s_to_T_restart**2 /= us%m_to_L_restart * us%s_to_T**2) ) then
1286  accel_rescale = (us%m_to_L * us%s_to_T_restart**2) / (us%m_to_L_restart * us%s_to_T**2)
1287  do k=1,nz ; do j=js,je ; do i=g%IscB,g%IecB
1288  cs%diffu(i,j,k) = accel_rescale * cs%diffu(i,j,k)
1289  enddo ; enddo ; enddo
1290  do k=1,nz ; do j=g%JscB,g%JecB ; do i=is,ie
1291  cs%diffv(i,j,k) = accel_rescale * cs%diffv(i,j,k)
1292  enddo ; enddo ; enddo
1293  endif
1294  endif
1295 
1296  if (.not. query_initialized(cs%u_av,"u2", restart_cs) .or. &
1297  .not. query_initialized(cs%u_av,"v2", restart_cs)) then
1298  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
1299  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
1300  elseif ( (us%s_to_T_restart * us%m_to_L_restart /= 0.0) .and. &
1301  ((us%m_to_L * us%s_to_T_restart) /= (us%m_to_L_restart * us%s_to_T)) ) then
1302  vel_rescale = (us%m_to_L * us%s_to_T_restart) / (us%m_to_L_restart * us%s_to_T)
1303  do k=1,nz ; do j=jsd,jed ; do i=isdb,iedb ; cs%u_av(i,j,k) = vel_rescale * cs%u_av(i,j,k) ; enddo ; enddo ; enddo
1304  do k=1,nz ; do j=jsdb,jedb ; do i=isd,ied ; cs%v_av(i,j,k) = vel_rescale * cs%v_av(i,j,k) ; enddo ; enddo ; enddo
1305  endif
1306 
1307  ! This call is just here to initialize uh and vh.
1308  if (.not. query_initialized(uh,"uh",restart_cs) .or. &
1309  .not. query_initialized(vh,"vh",restart_cs)) then
1310  do k=1,nz ; do j=jsd,jed ; do i=isd,ied ; h_tmp(i,j,k) = h(i,j,k) ; enddo ; enddo ; enddo
1311  call continuity(u, v, h, h_tmp, uh, vh, dt, g, gv, us, cs%continuity_CSp, obc=cs%OBC)
1312  call pass_var(h_tmp, g%Domain, clock=id_clock_pass_init)
1313  do k=1,nz ; do j=jsd,jed ; do i=isd,ied
1314  cs%h_av(i,j,k) = 0.5*(h(i,j,k) + h_tmp(i,j,k))
1315  enddo ; enddo ; enddo
1316  else
1317  if (.not. query_initialized(cs%h_av,"h2",restart_cs)) then
1318  cs%h_av(:,:,:) = h(:,:,:)
1319  elseif ((gv%m_to_H_restart /= 0.0) .and. (gv%m_to_H_restart /= gv%m_to_H)) then
1320  h_rescale = gv%m_to_H / gv%m_to_H_restart
1321  do k=1,nz ; do j=js,je ; do i=is,ie ; cs%h_av(i,j,k) = h_rescale * cs%h_av(i,j,k) ; enddo ; enddo ; enddo
1322  endif
1323  if ( (gv%m_to_H_restart * us%s_to_T_restart * us%m_to_L_restart /= 0.0) .and. &
1324  ((gv%m_to_H * us%m_to_L**2 * us%s_to_T_restart) /= &
1325  (gv%m_to_H_restart * us%m_to_L_restart**2 * us%s_to_T)) ) then
1326  uh_rescale = (gv%m_to_H * us%m_to_L**2 * us%s_to_T_restart) / &
1327  (gv%m_to_H_restart * us%m_to_L_restart**2 * us%s_to_T)
1328  do k=1,nz ; do j=js,je ; do i=g%IscB,g%IecB ; uh(i,j,k) = uh_rescale * uh(i,j,k) ; enddo ; enddo ; enddo
1329  do k=1,nz ; do j=g%JscB,g%JecB ; do i=is,ie ; vh(i,j,k) = uh_rescale * vh(i,j,k) ; enddo ; enddo ; enddo
1330  endif
1331  endif
1332 
1333  call cpu_clock_begin(id_clock_pass_init)
1334  call create_group_pass(pass_av_h_uvh, cs%u_av, cs%v_av, g%Domain, halo=2)
1335  call create_group_pass(pass_av_h_uvh, cs%h_av, g%Domain, halo=2)
1336  call create_group_pass(pass_av_h_uvh, uh, vh, g%Domain, halo=2)
1337  call do_group_pass(pass_av_h_uvh, g%Domain)
1338  call cpu_clock_end(id_clock_pass_init)
1339 
1340  flux_units = get_flux_units(gv)
1341  h_convert = gv%H_to_m ; if (.not.gv%Boussinesq) h_convert = gv%H_to_kg_m2
1342  cs%id_uh = register_diag_field('ocean_model', 'uh', diag%axesCuL, time, &
1343  'Zonal Thickness Flux', flux_units, y_cell_method='sum', v_extensive=.true., &
1344  conversion=h_convert*us%L_to_m**2*us%s_to_T)
1345  cs%id_vh = register_diag_field('ocean_model', 'vh', diag%axesCvL, time, &
1346  'Meridional Thickness Flux', flux_units, x_cell_method='sum', v_extensive=.true., &
1347  conversion=h_convert*us%L_to_m**2*us%s_to_T)
1348 
1349  cs%id_CAu = register_diag_field('ocean_model', 'CAu', diag%axesCuL, time, &
1350  'Zonal Coriolis and Advective Acceleration', 'm s-2', conversion=us%L_T2_to_m_s2)
1351  cs%id_CAv = register_diag_field('ocean_model', 'CAv', diag%axesCvL, time, &
1352  'Meridional Coriolis and Advective Acceleration', 'm s-2', conversion=us%L_T2_to_m_s2)
1353  cs%id_PFu = register_diag_field('ocean_model', 'PFu', diag%axesCuL, time, &
1354  'Zonal Pressure Force Acceleration', 'm s-2', conversion=us%L_T2_to_m_s2)
1355  cs%id_PFv = register_diag_field('ocean_model', 'PFv', diag%axesCvL, time, &
1356  'Meridional Pressure Force Acceleration', 'm s-2', conversion=us%L_T2_to_m_s2)
1357 
1358  !CS%id_hf_PFu = register_diag_field('ocean_model', 'hf_PFu', diag%axesCuL, Time, &
1359  ! 'Fractional Thickness-weighted Zonal Pressure Force Acceleration', 'm s-2', &
1360  ! v_extensive=.true., conversion=US%L_T2_to_m_s2)
1361  !if(CS%id_hf_PFu > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_u,IsdB,IedB,jsd,jed,nz)
1362 
1363  !CS%id_hf_PFv = register_diag_field('ocean_model', 'hf_PFv', diag%axesCvL, Time, &
1364  ! 'Fractional Thickness-weighted Meridional Pressure Force Acceleration', 'm s-2', &
1365  ! v_extensive=.true., conversion=US%L_T2_to_m_s2)
1366  !if(CS%id_hf_PFv > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
1367 
1368  !CS%id_hf_CAu = register_diag_field('ocean_model', 'hf_CAu', diag%axesCuL, Time, &
1369  ! 'Fractional Thickness-weighted Zonal Coriolis and Advective Acceleration', 'm s-2', &
1370  ! v_extensive=.true., conversion=US%L_T2_to_m_s2)
1371  !if(CS%id_hf_CAu > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_u,IsdB,IedB,jsd,jed,nz)
1372 
1373  !CS%id_hf_CAv = register_diag_field('ocean_model', 'hf_CAv', diag%axesCvL, Time, &
1374  ! 'Fractional Thickness-weighted Meridional Coriolis and Advective Acceleration', 'm s-2', &
1375  ! v_extensive=.true., conversion=US%L_T2_to_m_s2)
1376  !if(CS%id_hf_CAv > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
1377 
1378  cs%id_hf_PFu_2d = register_diag_field('ocean_model', 'hf_PFu_2d', diag%axesCu1, time, &
1379  'Depth-sum Fractional Thickness-weighted Zonal Pressure Force Acceleration', 'm s-2', &
1380  conversion=us%L_T2_to_m_s2)
1381  if(cs%id_hf_PFu_2d > 0) call safe_alloc_ptr(cs%ADp%diag_hfrac_u,isdb,iedb,jsd,jed,nz)
1382 
1383  cs%id_hf_PFv_2d = register_diag_field('ocean_model', 'hf_PFv_2d', diag%axesCv1, time, &
1384  'Depth-sum Fractional Thickness-weighted Meridional Pressure Force Acceleration', 'm s-2', &
1385  conversion=us%L_T2_to_m_s2)
1386  if(cs%id_hf_PFv_2d > 0) call safe_alloc_ptr(cs%ADp%diag_hfrac_v,isd,ied,jsd,jedb,nz)
1387 
1388  cs%id_hf_CAu_2d = register_diag_field('ocean_model', 'hf_CAu_2d', diag%axesCu1, time, &
1389  'Depth-sum Fractional Thickness-weighted Zonal Coriolis and Advective Acceleration', 'm s-2', &
1390  conversion=us%L_T2_to_m_s2)
1391  if(cs%id_hf_CAu_2d > 0) call safe_alloc_ptr(cs%ADp%diag_hfrac_u,isdb,iedb,jsd,jed,nz)
1392 
1393  cs%id_hf_CAv_2d = register_diag_field('ocean_model', 'hf_CAv_2d', diag%axesCv1, time, &
1394  'Depth-sum Fractional Thickness-weighted Meridional Coriolis and Advective Acceleration', 'm s-2', &
1395  conversion=us%L_T2_to_m_s2)
1396  if(cs%id_hf_CAv_2d > 0) call safe_alloc_ptr(cs%ADp%diag_hfrac_v,isd,ied,jsd,jedb,nz)
1397 
1398  cs%id_uav = register_diag_field('ocean_model', 'uav', diag%axesCuL, time, &
1399  'Barotropic-step Averaged Zonal Velocity', 'm s-1', conversion=us%L_T_to_m_s)
1400  cs%id_vav = register_diag_field('ocean_model', 'vav', diag%axesCvL, time, &
1401  'Barotropic-step Averaged Meridional Velocity', 'm s-1', conversion=us%L_T_to_m_s)
1402 
1403  cs%id_u_BT_accel = register_diag_field('ocean_model', 'u_BT_accel', diag%axesCuL, time, &
1404  'Barotropic Anomaly Zonal Acceleration', 'm s-2', conversion=us%L_T2_to_m_s2)
1405  cs%id_v_BT_accel = register_diag_field('ocean_model', 'v_BT_accel', diag%axesCvL, time, &
1406  'Barotropic Anomaly Meridional Acceleration', 'm s-2', conversion=us%L_T2_to_m_s2)
1407 
1408  !CS%id_hf_u_BT_accel = register_diag_field('ocean_model', 'hf_u_BT_accel', diag%axesCuL, Time, &
1409  ! 'Fractional Thickness-weighted Barotropic Anomaly Zonal Acceleration', 'm s-2', &
1410  ! v_extensive=.true., conversion=US%L_T2_to_m_s2)
1411  !if(CS%id_hf_u_BT_accel > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_u,IsdB,IedB,jsd,jed,nz)
1412 
1413  !CS%id_hf_v_BT_accel = register_diag_field('ocean_model', 'hf_v_BT_accel', diag%axesCvL, Time, &
1414  ! 'Fractional Thickness-weighted Barotropic Anomaly Meridional Acceleration', 'm s-2', &
1415  ! v_extensive=.true., conversion=US%L_T2_to_m_s2)
1416  !if(CS%id_hf_v_BT_accel > 0) call safe_alloc_ptr(CS%ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
1417 
1418  cs%id_hf_u_BT_accel_2d = register_diag_field('ocean_model', 'hf_u_BT_accel_2d', diag%axesCu1, time, &
1419  'Depth-sum Fractional Thickness-weighted Barotropic Anomaly Zonal Acceleration', 'm s-2', &
1420  conversion=us%L_T2_to_m_s2)
1421  if(cs%id_hf_u_BT_accel_2d > 0) call safe_alloc_ptr(cs%ADp%diag_hfrac_u,isdb,iedb,jsd,jed,nz)
1422 
1423  cs%id_hf_v_BT_accel_2d = register_diag_field('ocean_model', 'hf_v_BT_accel_2d', diag%axesCv1, time, &
1424  'Depth-sum Fractional Thickness-weighted Barotropic Anomaly Meridional Acceleration', 'm s-2', &
1425  conversion=us%L_T2_to_m_s2)
1426  if(cs%id_hf_v_BT_accel_2d > 0) call safe_alloc_ptr(cs%ADp%diag_hfrac_v,isd,ied,jsd,jedb,nz)
1427 
1428  id_clock_cor = cpu_clock_id('(Ocean Coriolis & mom advection)', grain=clock_module)
1429  id_clock_continuity = cpu_clock_id('(Ocean continuity equation)', grain=clock_module)
1430  id_clock_pres = cpu_clock_id('(Ocean pressure force)', grain=clock_module)
1431  id_clock_vertvisc = cpu_clock_id('(Ocean vertical viscosity)', grain=clock_module)
1432  id_clock_horvisc = cpu_clock_id('(Ocean horizontal viscosity)', grain=clock_module)
1433  id_clock_mom_update = cpu_clock_id('(Ocean momentum increments)', grain=clock_module)
1434  id_clock_pass = cpu_clock_id('(Ocean message passing)', grain=clock_module)
1435  id_clock_btcalc = cpu_clock_id('(Ocean barotropic mode calc)', grain=clock_module)
1436  id_clock_btstep = cpu_clock_id('(Ocean barotropic mode stepping)', grain=clock_module)
1437  id_clock_btforce = cpu_clock_id('(Ocean barotropic forcing calc)', grain=clock_module)
1438 
1439 end subroutine initialize_dyn_split_rk2
1440 
1441 
1442 !> Close the dyn_split_RK2 module
1443 subroutine end_dyn_split_rk2(CS)
1444  type(mom_dyn_split_rk2_cs), pointer :: cs !< module control structure
1445 
1446  dealloc_(cs%diffu) ; dealloc_(cs%diffv)
1447  dealloc_(cs%CAu) ; dealloc_(cs%CAv)
1448  dealloc_(cs%PFu) ; dealloc_(cs%PFv)
1449 
1450  if (associated(cs%taux_bot)) deallocate(cs%taux_bot)
1451  if (associated(cs%tauy_bot)) deallocate(cs%tauy_bot)
1452  dealloc_(cs%uhbt) ; dealloc_(cs%vhbt)
1453  dealloc_(cs%u_accel_bt) ; dealloc_(cs%v_accel_bt)
1454  dealloc_(cs%visc_rem_u) ; dealloc_(cs%visc_rem_v)
1455 
1456  dealloc_(cs%eta) ; dealloc_(cs%eta_PF) ; dealloc_(cs%pbce)
1457  dealloc_(cs%h_av) ; dealloc_(cs%u_av) ; dealloc_(cs%v_av)
1458 
1459  call dealloc_bt_cont_type(cs%BT_cont)
1460 
1461  deallocate(cs)
1462 end subroutine end_dyn_split_rk2
1463 
1464 
1465 !> \namespace mom_dynamics_split_rk2
1466 !!
1467 !! This file time steps the adiabatic dynamic core by splitting
1468 !! between baroclinic and barotropic modes. It uses a pseudo-second order
1469 !! Runge-Kutta time stepping scheme for the baroclinic momentum
1470 !! equation and a forward-backward coupling between the baroclinic
1471 !! momentum and continuity equations. This split time-stepping
1472 !! scheme is described in detail in Hallberg (JCP, 1997). Additional
1473 !! issues related to exact tracer conservation and how to
1474 !! ensure consistency between the barotropic and layered estimates
1475 !! of the free surface height are described in Hallberg and
1476 !! Adcroft (Ocean Modelling, 2009). This was the time stepping code
1477 !! that is used for most GOLD applications, including GFDL's ESM2G
1478 !! Earth system model, and all of the examples provided with the
1479 !! MOM code (although several of these solutions are routinely
1480 !! verified by comparison with the slower unsplit schemes).
1481 !!
1482 !! The subroutine step_MOM_dyn_split_RK2 actually does the time
1483 !! stepping, while register_restarts_dyn_split_RK2 sets the fields
1484 !! that are found in a full restart file with this scheme, and
1485 !! initialize_dyn_split_RK2 initializes the cpu clocks that are
1486 !! used in this module. For largely historical reasons, this module
1487 !! does not have its own control structure, but shares the same
1488 !! control structure with MOM.F90 and the other MOM_dynamics_...
1489 !! modules.
1490 
1491 end module mom_dynamics_split_rk2
mom_interface_heights::find_eta
Calculates the heights of sruface or all interfaces from layer thicknesses.
Definition: MOM_interface_heights.F90:21
mom_dynamics_split_rk2::mom_dyn_split_rk2_cs
MOM_dynamics_split_RK2 module control structure.
Definition: MOM_dynamics_split_RK2.F90:71
mom_time_manager
Wraps the FMS time manager functions.
Definition: MOM_time_manager.F90:2
mom_ale
This module contains the main regridding routines.
Definition: MOM_ALE.F90:9
mom_boundary_update::update_obc_cs
The control structure for the MOM_boundary_update module.
Definition: MOM_boundary_update.F90:37
mom_forcing_type
This module implements boundary forcing for MOM6.
Definition: MOM_forcing_type.F90:2
mom_open_boundary::ocean_obc_type
Open-boundary data.
Definition: MOM_open_boundary.F90:218
mom_hor_visc::hor_visc_cs
Control structure for horizontal viscosity.
Definition: MOM_hor_visc.F90:32
mom_interface_heights
Functions for calculating interface heights, including free surface height.
Definition: MOM_interface_heights.F90:2
mom_hor_visc
Calculates horizontal viscosity and viscous stresses.
Definition: MOM_hor_visc.F90:2
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:26
mom_diag_mediator::diag_ctrl
The following data type a list of diagnostic fields an their variants, as well as variables that cont...
Definition: MOM_diag_mediator.F90:240
mom_continuity
Solve the layer continuity equation.
Definition: MOM_continuity.F90:2
mom_boundary_update
Controls where open boundary conditions are applied.
Definition: MOM_boundary_update.F90:3
mom_get_input
Reads the only Fortran name list needed to boot-strap the model.
Definition: MOM_get_input.F90:6
mom_file_parser::param_file_type
A structure that can be parsed to read and document run-time parameters.
Definition: MOM_file_parser.F90:54
mom_vert_friction
Implements vertical viscosity (vertvisc)
Definition: MOM_vert_friction.F90:2
mom_grid
Provides the ocean grid type.
Definition: MOM_grid.F90:2
mom_variables::vertvisc_type
Vertical viscosities, drag coefficients, and related fields.
Definition: MOM_variables.F90:208
mom_cpu_clock
Wraps the MPP cpu clock functions.
Definition: MOM_cpu_clock.F90:2
mom_restart::register_restart_field
Register fields for restarts.
Definition: MOM_restart.F90:111
mom_coriolisadv
Accelerations due to the Coriolis force and momentum advection.
Definition: MOM_CoriolisAdv.F90:2
mom_set_visc::set_visc_cs
Control structure for MOM_set_visc.
Definition: MOM_set_viscosity.F90:45
mom_pressureforce
A thin wrapper for Boussinesq/non-Boussinesq forms of the pressure force calculation.
Definition: MOM_PressureForce.F90:2
mom_meke_types::meke_type
This type is used to exchange information related to the MEKE calculations.
Definition: MOM_MEKE_types.F90:8
mom_io
This module contains I/O framework code.
Definition: MOM_io.F90:2
mom_set_visc
Calculates various values related to the bottom boundary layer, such as the viscosity and thickness o...
Definition: MOM_set_viscosity.F90:3
mom_file_parser
The MOM6 facility to parse input files for runtime parameters.
Definition: MOM_file_parser.F90:2
mom_hor_index
Defines the horizontal index type (hor_index_type) used for providing index ranges.
Definition: MOM_hor_index.F90:2
mom_variables::accel_diag_ptrs
Pointers to arrays with accelerations, which can later be used for derived diagnostics,...
Definition: MOM_variables.F90:161
mom_restart::register_restart_pair
Register a pair of restart fieilds whose rotations map onto each other.
Definition: MOM_restart.F90:120
mom_lateral_mixing_coeffs
Variable mixing coefficients.
Definition: MOM_lateral_mixing_coeffs.F90:2
mom_hor_index::hor_index_type
Container for horizontal index ranges for data, computational and global domains.
Definition: MOM_hor_index.F90:16
mom_coriolisadv::coriolisadv_cs
Control structure for mom_coriolisadv.
Definition: MOM_CoriolisAdv.F90:27
mom_forcing_type::mech_forcing
Structure that contains pointers to the mechanical forcing at the surface used to drive the liquid oc...
Definition: MOM_forcing_type.F90:204
mom_lateral_mixing_coeffs::varmix_cs
Variable mixing coefficients.
Definition: MOM_lateral_mixing_coeffs.F90:27
mom_barotropic
Baropotric solver.
Definition: MOM_barotropic.F90:2
mom_wave_interface
Interface for surface waves.
Definition: MOM_wave_interface.F90:2
mom_diag_mediator::post_data
Make a diagnostic available for averaging or output.
Definition: MOM_diag_mediator.F90:70
mom_restart::mom_restart_cs
A restart registry and the control structure for restarts.
Definition: MOM_restart.F90:75
mom_variables::bt_cont_type
Container for information about the summed layer transports and how they will vary as the barotropic ...
Definition: MOM_variables.F90:263
mom_unit_scaling::unit_scale_type
Describes various unit conversion factors.
Definition: MOM_unit_scaling.F90:14
mom_continuity::continuity_cs
Control structure for mom_continuity.
Definition: MOM_continuity.F90:26
mom_checksum_packages
Provides routines that do checksums of groups of MOM variables.
Definition: MOM_checksum_packages.F90:2
mom_tidal_forcing
Tidal contributions to geopotential.
Definition: MOM_tidal_forcing.F90:2
mom_unit_scaling
Provides a transparent unit rescaling type to facilitate dimensional consistency testing.
Definition: MOM_unit_scaling.F90:2
mom_domains
Describes the decomposed MOM domain and has routines for communications across PEs.
Definition: MOM_domains.F90:2
mom_thickness_diffuse
Thickness diffusion (or Gent McWilliams)
Definition: MOM_thickness_diffuse.F90:2
mom_diag_mediator
The subroutines here provide convenient wrappers to the fms diag_manager interfaces with additional d...
Definition: MOM_diag_mediator.F90:3
mom_variables::cont_diag_ptrs
Pointers to arrays with transports, which can later be used for derived diagnostics,...
Definition: MOM_variables.F90:193
mom_error_handler
Routines for error handling and I/O management.
Definition: MOM_error_handler.F90:2
mom_restart
The MOM6 facility for reading and writing restart files, and querying what has been read.
Definition: MOM_restart.F90:2
mom_variables::ocean_internal_state
Pointers to all of the prognostic variables allocated in MOM_variables.F90 and MOM....
Definition: MOM_variables.F90:132
mom_wave_interface::wave_parameters_cs
Container for all surface wave related parameters.
Definition: MOM_wave_interface.F90:47
mom_vert_friction::vertvisc_cs
The control structure with parameters and memory for the MOM_vert_friction module.
Definition: MOM_vert_friction.F90:39
mom_io::vardesc
Type for describing a variable, typically a tracer.
Definition: MOM_io.F90:53
mom_file_parser::log_version
An overloaded interface to log version information about modules.
Definition: MOM_file_parser.F90:109
mom_verticalgrid::verticalgrid_type
Describes the vertical ocean grid, including unit conversion factors.
Definition: MOM_verticalGrid.F90:24
mom_domains::create_group_pass
Set up a group of halo updates.
Definition: MOM_domains.F90:84
mom_debugging::check_redundant
Check for consistency between the duplicated points of a C-grid vector.
Definition: MOM_debugging.F90:33
mom_get_input::directories
Container for paths and parameter file names.
Definition: MOM_get_input.F90:20
mom_variables::thermo_var_ptrs
Pointers to an assortment of thermodynamic fields that may be available, including potential temperat...
Definition: MOM_variables.F90:80
mom_tidal_forcing::tidal_forcing_cs
The control structure for the MOM_tidal_forcing module.
Definition: MOM_tidal_forcing.F90:36
mom_thickness_diffuse::thickness_diffuse_cs
Control structure for thickness diffusion.
Definition: MOM_thickness_diffuse.F90:37
mom_restart::query_initialized
Indicate whether a field has been read from a restart file.
Definition: MOM_restart.F90:127
mom_barotropic::barotropic_cs
The barotropic stepping control stucture.
Definition: MOM_barotropic.F90:102
mom_open_boundary
Controls where open boundary conditions are applied.
Definition: MOM_open_boundary.F90:2
mom_verticalgrid
Provides a transparent vertical ocean grid type and supporting routines.
Definition: MOM_verticalGrid.F90:2
mom_ale::ale_cs
ALE control structure.
Definition: MOM_ALE.F90:63
mom_variables
Provides transparent structures with groups of MOM6 variables and supporting routines.
Definition: MOM_variables.F90:2
mom_pressureforce::pressureforce_cs
Pressure force control structure.
Definition: MOM_PressureForce.F90:28
mom_domains::pass_var
Do a halo update on an array.
Definition: MOM_domains.F90:54
mom_dynamics_split_rk2
Time step the adiabatic dynamic core of MOM using RK2 method.
Definition: MOM_dynamics_split_RK2.F90:2
mom_checksum_packages::mom_state_chksum
Write out checksums of the MOM6 state variables.
Definition: MOM_checksum_packages.F90:23
mom_file_parser::get_param
An overloaded interface to read and log the values of various types of parameters.
Definition: MOM_file_parser.F90:102
mom_debugging
Provides checksumming functions for debugging.
Definition: MOM_debugging.F90:7