MOM_dynamics_unsplit_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 steps the ocean dynamics with an unsplit quasi 2nd order Runge-Kutta scheme
7
8!********+*********+*********+*********+*********+*********+*********+**
9!* *
10!* By Alistair Adcroft and Robert Hallberg, 2010-2012 *
11!* *
12!* This file contains code that does the time-stepping of the *
13!* adiabatic dynamic core, in this case with a pseudo-second order *
14!* Runge-Kutta time stepping scheme for the momentum and a forward- *
15!* backward coupling between the momentum and continuity equations, *
16!* but without any splitting between the baroclinic and barotropic *
17!* modes. Apart from the lack of splitting, this is closely analogous *
18!* to the split time stepping scheme, and efforts have been taken to *
19!* ensure that for certain configurations (e.g., very short *
20!* baroclinic time steps, a single barotropic step per baroclinic *
21!* step, and particular choices about how to coupled the baroclinic *
22!* and barotropic solves, the two solutions reproduce each other. *
23!* Although this time stepping scheme is not very efficient with a *
24!* large number of layers, it is valuable for verifying the proper *
25!* behavior of the more complicated split time stepping scheme, and *
26!* is not too inefficient for use with only a few layers. *
27!* *
28!* The subroutine step_MOM_dyn_unsplit_RK2 actually does the time *
29!* stepping, while register_restarts_dyn_unsplit_RK2 sets the fields *
30!* that are found in a full restart file with this scheme, and *
31!* initialize_dyn_unsplit_RK2 initializes the cpu clocks that are *
32!* used in this module. For largely historical reasons, this module *
33!* does not have its own control structure, but shares the same *
34!* control structure with MOM.F90 and the other MOM_dynamics_... *
35!* modules. *
36!* *
37!* Macros written all in capital letters are defined in MOM_memory.h. *
38!* *
39!* A small fragment of the grid is shown below: *
40!* *
41!* j+1 x ^ x ^ x At x: q, CoriolisBu *
42!* j+1 > o > o > At ^: v, PFv, CAv, vh, diffv, tauy, vbt, vhtr *
43!* j x ^ x ^ x At >: u, PFu, CAu, uh, diffu, taux, ubt, uhtr *
44!* j > o > o > At o: h, bathyT, eta, T, S, tr *
45!* j-1 x ^ x ^ x *
46!* i-1 i i+1 *
47!* i i+1 *
48!* *
49!* The boundaries always run through q grid points (x). *
50!* *
51!********+*********+*********+*********+*********+*********+*********+**
52
56use mom_checksum_packages, only : mom_thermo_chksum, mom_state_chksum, mom_accel_chksum
57use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
58use mom_cpu_clock, only : clock_component, clock_subcomponent
59use mom_cpu_clock, only : clock_module_driver, clock_module, clock_routine
61use mom_diag_mediator, only : disable_averaging, post_data, safe_alloc_ptr
64use mom_domains, only : pass_var, pass_var_start, pass_var_complete
65use mom_domains, only : pass_vector, pass_vector_start, pass_vector_complete
66use mom_domains, only : to_south, to_west, to_all, cgrid_ne, scalar_pair
67use mom_error_handler, only : mom_error, mom_mesg, fatal, warning, is_root_pe
69use mom_file_parser, only : get_param, log_param, log_version, param_file_type
70use mom_get_input, only : directories
71use mom_time_manager, only : time_type, operator(+)
72use mom_time_manager, only : operator(-), operator(>), operator(*), operator(/)
73
74use mom_ale, only : ale_cs
77use mom_continuity, only : continuity, continuity_init, continuity_cs, continuity_stencil
79use mom_debugging, only : check_redundant
80use mom_grid, only : ocean_grid_type
83use mom_interface_heights, only : thickness_to_dz
85use mom_meke_types, only : meke_type
98
99implicit none ; private
100
101#include <MOM_memory.h>
102
103!> MOM_dynamics_unsplit_RK2 module control structure
104type, public :: mom_dyn_unsplit_rk2_cs ; private
105 real allocable_, dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: &
106 cau, & !< CAu = f*v - u.grad(u) [L T-2 ~> m s-2].
107 pfu, & !< PFu = -dM/dx [L T-2 ~> m s-2].
108 diffu !< Zonal acceleration due to convergence of the along-isopycnal stress tensor [L T-2 ~> m s-2].
109
110 real allocable_, dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: &
111 cav, & !< CAv = -f*u - u.grad(v) [L T-2 ~> m s-2].
112 pfv, & !< PFv = -dM/dy [L T-2 ~> m s-2].
113 diffv !< Meridional acceleration due to convergence of the along-isopycnal stress tensor [L T-2 ~> m s-2].
114
115 real, pointer, dimension(:,:) :: taux_bot => null() !< frictional x-bottom stress from the ocean
116 !! to the seafloor [R L Z T-2 ~> Pa]
117 real, pointer, dimension(:,:) :: tauy_bot => null() !< frictional y-bottom stress from the ocean
118 !! to the seafloor [R L Z T-2 ~> Pa]
119
120 real :: be !< A nondimensional number from 0.5 to 1 that controls
121 !! the backward weighting of the time stepping scheme [nondim].
122 real :: begw !< A nondimensional number from 0 to 1 that controls
123 !! the extent to which the treatment of gravity waves
124 !! is forward-backward (0) or simulated backward
125 !! Euler (1) [nondim]. 0 is often used.
126 logical :: dt_visc_bug !< If false, use the correct timestep in the calculation of the
127 !! turbulent mixed layer properties for viscosity. Otherwise if
128 !! this is true, an older incorrect setting is used.
129 logical :: debug !< If true, write verbose checksums for debugging purposes.
130 logical :: calculate_sal !< If true, calculate self-attraction and loading.
131 logical :: use_tides !< If true, tidal forcing is enabled.
132
133 logical :: module_is_initialized = .false. !< Record whether this module has been initialized.
134
135 !>@{ Diagnostic IDs
136 integer :: id_uh = -1, id_vh = -1
137 integer :: id_ueffa = -1, id_veffa = -1
138 integer :: id_pfu = -1, id_pfv = -1, id_cau = -1, id_cav = -1
139 !>@}
140
141 type(diag_ctrl), pointer :: diag => null() !< A structure that is used to
142 !! regulate the timing of diagnostic output.
143 type(accel_diag_ptrs), pointer :: adp => null() !< A structure pointing to the
144 !! accelerations in the momentum equations,
145 !! which can later be used to calculate
146 !! derived diagnostics like energy budgets.
147 type(cont_diag_ptrs), pointer :: cdp => null() !< A structure with pointers to
148 !! various terms in the continuity equations,
149 !! which can later be used to calculate
150 !! derived diagnostics like energy budgets.
151
152 ! The remainder of the structure points to child subroutines' control structures.
153 !> A pointer to the horizontal viscosity control structure
154 type(hor_visc_cs) :: hor_visc
155 !> A pointer to the continuity control structure
156 type(continuity_cs) :: continuity_csp
157 !> The CoriolisAdv control structure
158 type(coriolisadv_cs) :: coriolisadv
159 !> A pointer to the PressureForce control structure
160 type(pressureforce_cs) :: pressureforce_csp
161 !> A pointer to the vertvisc control structure
162 type(vertvisc_cs), pointer :: vertvisc_csp => null()
163 !> A pointer to the set_visc control structure
164 type(set_visc_cs), pointer :: set_visc_csp => null()
165 !> A pointer to the SAL control structure
166 type(sal_cs) :: sal_csp
167 !> A pointer to the tidal forcing control structure
168 type(tidal_forcing_cs) :: tides_csp
169 !> A pointer to the ALE control structure.
170 type(ale_cs), pointer :: ale_csp => null()
171
172 type(ocean_obc_type), pointer :: obc => null() !< A pointer to an open boundary
173 !! condition type that specifies whether, where, and what open boundary
174 !! conditions are used. If no open BCs are used, this pointer stays
175 !! nullified. Flather OBCs use open boundary_CS as well.
176 !> A pointer to the update_OBC control structure
177 type(update_obc_cs), pointer :: update_obc_csp => null()
178
180
181
184
185!>@{ CPU time clock IDs
186integer :: id_clock_cor, id_clock_pres, id_clock_vertvisc
187integer :: id_clock_horvisc, id_clock_continuity, id_clock_mom_update
188integer :: id_clock_pass, id_clock_pass_init
189!>@}
190
191contains
192
193! =============================================================================
194
195!> Step the MOM6 dynamics using an unsplit quasi-2nd order Runge-Kutta scheme
196subroutine step_mom_dyn_unsplit_rk2(u_in, v_in, h_in, tv, visc, Time_local, dt, forces, &
197 p_surf_begin, p_surf_end, uh, vh, uhtr, vhtr, eta_av, G, GV, US, CS, &
198 VarMix, MEKE, pbv, STOCH)
199 type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure.
200 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
201 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
202 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: u_in !< The input and output zonal
203 !! velocity [L T-1 ~> m s-1].
204 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(inout) :: v_in !< The input and output meridional
205 !! velocity [L T-1 ~> m s-1].
206 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h_in !< The input and output layer thicknesses,
207 !! [H ~> m or kg m-2], depending on whether
208 !! the Boussinesq approximation is made.
209 type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various
210 !! thermodynamic variables.
211 type(vertvisc_type), intent(inout) :: visc !< A structure containing vertical
212 !! viscosities, bottom drag
213 !! viscosities, and related fields.
214 type(time_type), intent(in) :: time_local !< The model time at the end of
215 !! the time step.
216 real, intent(in) :: dt !< The baroclinic dynamics time step [T ~> s].
217 type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces
218 real, dimension(:,:), pointer :: p_surf_begin !< A pointer (perhaps NULL) to
219 !! the surface pressure at the beginning
220 !! of this dynamic step [R L2 T-2 ~> Pa].
221 real, dimension(:,:), pointer :: p_surf_end !< A pointer (perhaps NULL) to
222 !! the surface pressure at the end of
223 !! this dynamic step [R L2 T-2 ~> Pa].
224 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: uh !< The zonal volume or mass transport
225 !! [H L2 T-1 ~> m3 s-1 or kg s-1].
226 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(inout) :: vh !< The meridional volume or mass
227 !! transport [H L2 T-1 ~> m3 s-1 or kg s-1].
228 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: uhtr !< The accumulated zonal volume or
229 !! mass transport since the last
230 !! tracer advection [H L2 ~> m3 or kg].
231 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(inout) :: vhtr !< The accumulated meridional volume
232 !! or mass transport since the last
233 !! tracer advection [H L2 ~> m3 or kg].
234 real, dimension(SZI_(G),SZJ_(G)), intent(out) :: eta_av !< The time-mean free surface height
235 !! or column mass [H ~> m or kg m-2].
236 type(mom_dyn_unsplit_rk2_cs), pointer :: cs !< The control structure set up by
237 !! initialize_dyn_unsplit_RK2.
238 type(varmix_cs), intent(inout) :: varmix !< Variable mixing control structure
239 type(meke_type), intent(inout) :: meke !< MEKE fields
240 !! fields related to the Mesoscale
241 !! Eddy Kinetic Energy.
242 type(porous_barrier_type), intent(in) :: pbv !< porous barrier fractional cell metrics
243 type(stochastic_cs), intent(inout) :: stoch !< Stochastic control structure
244 ! Local variables
245 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_av ! Averaged layer thicknesses [H ~> m or kg m-2]
246 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: hp ! Predicted layer thicknesses [H ~> m or kg m-2]
247 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: dz ! Distance between the interfaces around a layer [Z ~> m]
248 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: up ! Predicted zonal velocities [L T-1 ~> m s-1]
249 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: vp ! Predicted meridional velocities [L T-1 ~> m s-1]
250 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: ueffa ! Effective Area of U-Faces [H L ~> m2]
251 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: veffa ! Effective Area of V-Faces [H L ~> m2]
252 real, dimension(:,:), pointer :: p_surf => null() ! A pointer to the surface pressure [R L2 T-2 ~> Pa]
253 real :: dt_pred ! The time step for the predictor part of the baroclinic time stepping [T ~> s]
254 real :: dt_visc ! The time step for a part of the update due to viscosity [T ~> s]
255 logical :: dyn_p_surf
256 integer :: i, j, k, is, ie, js, je, isq, ieq, jsq, jeq, nz
257 integer :: cor_stencil ! Stencil size for Coriolis schemes [nondim]
258 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
259 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
260 dt_pred = dt * cs%BE
261 cor_stencil = coriolisadv_stencil(cs%CoriolisAdv)
262
263 h_av(:,:,:) = 0 ; hp(:,:,:) = 0
264 up(:,:,:) = 0
265 vp(:,:,:) = 0
266
267 dyn_p_surf = associated(p_surf_begin) .and. associated(p_surf_end)
268 if (dyn_p_surf) then
269 call safe_alloc_ptr(p_surf,g%isd,g%ied,g%jsd,g%jed) ; p_surf(:,:) = 0.0
270 else
271 p_surf => forces%p_surf
272 endif
273
274! Runge-Kutta second order accurate two step scheme is used to step
275! all of the fields except h. h is stepped separately.
276
277 if (cs%debug) then
278 call mom_state_chksum("Start Predictor ", u_in, v_in, h_in, uh, vh, g, gv, us)
279 endif
280
281! diffu = horizontal viscosity terms (u,h)
282 call enable_averages(dt,time_local, cs%diag)
283 call cpu_clock_begin(id_clock_horvisc)
284 call horizontal_viscosity(u_in, v_in, h_in, uh, vh, cs%diffu, cs%diffv, meke, varmix, &
285 g, gv, us, cs%hor_visc, tv, dt, stoch=stoch)
286 call cpu_clock_end(id_clock_horvisc)
287 call disable_averaging(cs%diag)
288 call pass_vector(cs%diffu, cs%diffv, g%Domain, clock=id_clock_pass)
289
290! This continuity step is solely for the Coroilis terms, specifically in the
291! denominator of PV and in the mass transport or PV.
292! uh = u[n-1]*h[n-1/2]
293! hp = h[n-1/2] + dt/2 div . uh
294 call cpu_clock_begin(id_clock_continuity)
295 ! This is a duplicate calculation of the last continuity from the previous step
296 ! and could/should be optimized out. -AJA
297 call continuity(u_in, v_in, h_in, hp, uh, vh, dt_pred, g, gv, us, cs%continuity_CSp, cs%OBC, pbv)
298 call cpu_clock_end(id_clock_continuity)
299 call pass_var(hp, g%Domain, halo=cor_stencil, clock=id_clock_pass)
300 call pass_vector(uh, vh, g%Domain, halo=cor_stencil, clock=id_clock_pass)
301 if (cor_stencil > 2) then
302 call pass_vector(u_in, v_in, g%Domain, halo=cor_stencil, clock=id_clock_pass)
303 endif
304
305! h_av = (h + hp)/2 (used in PV denominator)
306 call cpu_clock_begin(id_clock_mom_update)
307 do k=1,nz ; do j=js-cor_stencil,je+cor_stencil ; do i=is-cor_stencil,ie+cor_stencil
308 h_av(i,j,k) = (h_in(i,j,k) + hp(i,j,k)) * 0.5
309 enddo ; enddo ; enddo
310 call cpu_clock_end(id_clock_mom_update)
311
312! CAu = -(f+zeta)/h_av vh + d/dx KE (function of u[n-1] and uh[n-1])
313 call cpu_clock_begin(id_clock_cor)
314 call coradcalc(u_in, v_in, h_av, uh, vh, cs%CAu, cs%CAv, cs%OBC, cs%ADp, &
315 g, gv, us, cs%CoriolisAdv, pbv)
316 call cpu_clock_end(id_clock_cor)
317
318! PFu = d/dx M(h_av,T,S) (function of h[n-1/2])
319 call cpu_clock_begin(id_clock_pres)
320 if (dyn_p_surf) then ; do j=js-2,je+2 ; do i=is-2,ie+2
321 p_surf(i,j) = 0.5*p_surf_begin(i,j) + 0.5*p_surf_end(i,j)
322 enddo ; enddo ; endif
323 call pressureforce(h_in, tv, cs%PFu, cs%PFv, g, gv, us, cs%PressureForce_CSp, cs%ALE_CSp, cs%ADp, p_surf)
324 call cpu_clock_end(id_clock_pres)
325 call pass_vector(cs%PFu, cs%PFv, g%Domain, clock=id_clock_pass)
326 call pass_vector(cs%CAu, cs%CAv, g%Domain, clock=id_clock_pass)
327
328 if (associated(cs%OBC)) then ; if (cs%OBC%update_OBC) then
329 call update_obc_data(cs%OBC, g, gv, us, tv, h_in, cs%update_OBC_CSp, time_local)
330 endif ; endif
331 if (associated(cs%OBC)) then
332 call open_boundary_zero_normal_flow(cs%OBC, g, gv, cs%PFu, cs%PFv)
333 call open_boundary_zero_normal_flow(cs%OBC, g, gv, cs%CAu, cs%CAv)
334 call open_boundary_zero_normal_flow(cs%OBC, g, gv, cs%diffu, cs%diffv)
335 endif
336
337! up+[n-1/2] = u[n-1] + dt_pred * (PFu + CAu)
338 call cpu_clock_begin(id_clock_mom_update)
339 do k=1,nz ; do j=js,je ; do i=isq,ieq
340 up(i,j,k) = g%mask2dCu(i,j) * (u_in(i,j,k) + dt_pred * &
341 ((cs%PFu(i,j,k) + cs%CAu(i,j,k)) + cs%diffu(i,j,k)))
342 enddo ; enddo ; enddo
343 do k=1,nz ; do j=jsq,jeq ; do i=is,ie
344 vp(i,j,k) = g%mask2dCv(i,j) * (v_in(i,j,k) + dt_pred * &
345 ((cs%PFv(i,j,k) + cs%CAv(i,j,k)) + cs%diffv(i,j,k)))
346 enddo ; enddo ; enddo
347 call cpu_clock_end(id_clock_mom_update)
348
349 if (cs%debug) &
350 call mom_accel_chksum("Predictor 1 accel", cs%CAu, cs%CAv, cs%PFu, cs%PFv,&
351 cs%diffu, cs%diffv, g, gv, us)
352
353 ! up[n-1/2] <- up*[n-1/2] + dt/2 d/dz visc d/dz up[n-1/2]
354 call cpu_clock_begin(id_clock_vertvisc)
355 call enable_averages(dt, time_local, cs%diag)
356 dt_visc = dt ; if (cs%dt_visc_bug) dt_visc = dt_pred
357 call set_viscous_ml(u_in, v_in, h_av, tv, forces, visc, dt_visc, g, gv, us, cs%set_visc_CSp)
358 call disable_averaging(cs%diag)
359
360 call thickness_to_dz(h_av, tv, dz, g, gv, us, halo_size=1)
361 call vertvisc_coef(up, vp, h_av, dz, forces, visc, tv, dt_pred, g, gv, us, cs%vertvisc_CSp, cs%OBC, varmix)
362 call vertvisc(up, vp, h_av, forces, visc, dt_pred, cs%OBC, cs%ADp, cs%CDp, &
363 g, gv, us, cs%vertvisc_CSp)
364 call cpu_clock_end(id_clock_vertvisc)
365 call pass_vector(up, vp, g%Domain, clock=id_clock_pass)
366
367! uh = up[n-1/2] * h[n-1/2]
368! h_av = h + dt div . uh
369 call cpu_clock_begin(id_clock_continuity)
370 call continuity(up, vp, h_in, hp, uh, vh, dt, g, gv, us, cs%continuity_CSp, cs%OBC, pbv)
371 call cpu_clock_end(id_clock_continuity)
372 call pass_var(hp, g%Domain, clock=id_clock_pass)
373 call pass_vector(uh, vh, g%Domain, clock=id_clock_pass)
374
375! h_av <- (h + hp)/2 (centered at n-1/2)
376 do k=1,nz ; do j=js-cor_stencil,je+cor_stencil ; do i=is-cor_stencil,ie+cor_stencil
377 h_av(i,j,k) = (h_in(i,j,k) + hp(i,j,k)) * 0.5
378 enddo ; enddo ; enddo
379
380 if (cs%debug) &
381 call mom_state_chksum("Predictor 1", up, vp, h_av, uh, vh, g, gv, us)
382
383! CAu = -(f+zeta(up))/h_av vh + d/dx KE(up) (function of up[n-1/2], h[n-1/2])
384 call cpu_clock_begin(id_clock_cor)
385 call coradcalc(up, vp, h_av, uh, vh, cs%CAu, cs%CAv, cs%OBC, cs%ADp, &
386 g, gv, us, cs%CoriolisAdv, pbv)
387 call cpu_clock_end(id_clock_cor)
388 if (associated(cs%OBC)) then
389 call open_boundary_zero_normal_flow(cs%OBC, g, gv, cs%CAu, cs%CAv)
390 endif
391
392! call enable_averages(dt, Time_local, CS%diag) ?????????????????????/
393
394! up* = u[n] + (1+gamma) * dt * ( PFu + CAu ) Extrapolated for damping
395! u*[n+1] = u[n] + dt * ( PFu + CAu )
396 do k=1,nz ; do j=js,je ; do i=isq,ieq
397 up(i,j,k) = g%mask2dCu(i,j) * (u_in(i,j,k) + dt * (1.+cs%begw) * &
398 ((cs%PFu(i,j,k) + cs%CAu(i,j,k)) + cs%diffu(i,j,k)))
399 u_in(i,j,k) = g%mask2dCu(i,j) * (u_in(i,j,k) + dt * &
400 ((cs%PFu(i,j,k) + cs%CAu(i,j,k)) + cs%diffu(i,j,k)))
401 enddo ; enddo ; enddo
402 do k=1,nz ; do j=jsq,jeq ; do i=is,ie
403 vp(i,j,k) = g%mask2dCv(i,j) * (v_in(i,j,k) + dt * (1.+cs%begw) * &
404 ((cs%PFv(i,j,k) + cs%CAv(i,j,k)) + cs%diffv(i,j,k)))
405 v_in(i,j,k) = g%mask2dCv(i,j) * (v_in(i,j,k) + dt * &
406 ((cs%PFv(i,j,k) + cs%CAv(i,j,k)) + cs%diffv(i,j,k)))
407 enddo ; enddo ; enddo
408
409! up[n] <- up* + dt d/dz visc d/dz up
410! u[n] <- u*[n] + dt d/dz visc d/dz u[n]
411 call cpu_clock_begin(id_clock_vertvisc)
412 call thickness_to_dz(h_av, tv, dz, g, gv, us, halo_size=1)
413 call vertvisc_coef(up, vp, h_av, dz, forces, visc, tv, dt, g, gv, us, cs%vertvisc_CSp, cs%OBC, varmix)
414 call vertvisc(up, vp, h_av, forces, visc, dt, cs%OBC, cs%ADp, cs%CDp, &
415 g, gv, us, cs%vertvisc_CSp, cs%taux_bot, cs%tauy_bot)
416 call vertvisc_coef(u_in, v_in, h_av, dz, forces, visc, tv, dt, g, gv, us, cs%vertvisc_CSp, cs%OBC, varmix)
417 call vertvisc(u_in, v_in, h_av, forces, visc, dt, cs%OBC, cs%ADp, cs%CDp,&
418 g, gv, us, cs%vertvisc_CSp, cs%taux_bot, cs%tauy_bot)
419 call cpu_clock_end(id_clock_vertvisc)
420 call pass_vector(up, vp, g%Domain, clock=id_clock_pass)
421 call pass_vector(u_in, v_in, g%Domain, clock=id_clock_pass)
422
423! uh = up[n] * h[n] (up[n] might be extrapolated to damp GWs)
424! h[n+1] = h[n] + dt div . uh
425 call cpu_clock_begin(id_clock_continuity)
426 call continuity(up, vp, h_in, h_in, uh, vh, dt, g, gv, us, cs%continuity_CSp, cs%OBC, pbv)
427 call cpu_clock_end(id_clock_continuity)
428 call pass_var(h_in, g%Domain, clock=id_clock_pass)
429 call pass_vector(uh, vh, g%Domain, clock=id_clock_pass)
430
431! Accumulate mass flux for tracer transport
432 do k=1,nz
433 do j=js-2,je+2 ; do i=isq-2,ieq+2
434 uhtr(i,j,k) = uhtr(i,j,k) + dt*uh(i,j,k)
435 enddo ; enddo
436 do j=jsq-2,jeq+2 ; do i=is-2,ie+2
437 vhtr(i,j,k) = vhtr(i,j,k) + dt*vh(i,j,k)
438 enddo ; enddo
439 enddo
440
441 if (cs%debug) then
442 call mom_state_chksum("Corrector", u_in, v_in, h_in, uh, vh, g, gv, us)
443 call mom_accel_chksum("Corrector accel", cs%CAu, cs%CAv, cs%PFu, cs%PFv, &
444 cs%diffu, cs%diffv, g, gv, us)
445 endif
446
447 if (gv%Boussinesq) then
448 do j=js,je ; do i=is,ie ; eta_av(i,j) = -gv%Z_to_H*g%bathyT(i,j) ; enddo ; enddo
449 else
450 do j=js,je ; do i=is,ie ; eta_av(i,j) = 0.0 ; enddo ; enddo
451 endif
452 do k=1,nz ; do j=js,je ; do i=is,ie
453 eta_av(i,j) = eta_av(i,j) + h_av(i,j,k)
454 enddo ; enddo ; enddo
455
456 if (dyn_p_surf) deallocate(p_surf)
457
458! Here various terms used in to update the momentum equations are
459! offered for averaging.
460 if (cs%id_PFu > 0) call post_data(cs%id_PFu, cs%PFu, cs%diag)
461 if (cs%id_PFv > 0) call post_data(cs%id_PFv, cs%PFv, cs%diag)
462 if (cs%id_CAu > 0) call post_data(cs%id_CAu, cs%CAu, cs%diag)
463 if (cs%id_CAv > 0) call post_data(cs%id_CAv, cs%CAv, cs%diag)
464
465! Here the thickness fluxes are offered for averaging.
466 if (cs%id_uh > 0) call post_data(cs%id_uh, uh, cs%diag)
467 if (cs%id_vh > 0) call post_data(cs%id_vh, vh, cs%diag)
468
469! Calculate effective areas and post data
470 if (cs%id_ueffA > 0) then
471 ueffa(:,:,:) = 0
472 do k=1,nz ; do j=js,je ; do i=isq,ieq
473 if (abs(up(i,j,k)) > 0.) ueffa(i,j,k) = uh(i,j,k)/up(i,j,k)
474 enddo ; enddo ; enddo
475 call post_data(cs%id_ueffA, ueffa, cs%diag)
476 endif
477
478 if (cs%id_veffA > 0) then
479 veffa(:,:,:) = 0
480 do k=1,nz ; do j=jsq,jeq ; do i=is,ie
481 if (abs(vp(i,j,k)) > 0.) veffa(i,j,k) = vh(i,j,k)/vp(i,j,k)
482 enddo ; enddo ; enddo
483 call post_data(cs%id_veffA, veffa, cs%diag)
484 endif
485
486
487end subroutine step_mom_dyn_unsplit_rk2
488
489! =============================================================================
490
491!> Allocate the control structure for this module, allocates memory in it, and registers
492!! any auxiliary restart variables that are specific to the unsplit RK2 time stepping scheme.
493!!
494!! All variables registered here should have the ability to be recreated if they are not present
495!! in a restart file.
496subroutine register_restarts_dyn_unsplit_rk2(HI, GV, param_file, CS)
497 type(hor_index_type), intent(in) :: hi !< A horizontal index type structure.
498 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
499 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time
500 !! parameters.
501 type(mom_dyn_unsplit_rk2_cs), pointer :: cs !< The control structure set up by
502 !! initialize_dyn_unsplit_RK2.
503! This subroutine sets up any auxiliary restart variables that are specific
504! to the unsplit time stepping scheme. All variables registered here should
505! have the ability to be recreated if they are not present in a restart file.
506
507 ! Local variables
508 character(len=48) :: thickness_units, flux_units
509 integer :: isd, ied, jsd, jed, nz, isdb, iedb, jsdb, jedb
510 isd = hi%isd ; ied = hi%ied ; jsd = hi%jsd ; jed = hi%jed ; nz = gv%ke
511 isdb = hi%IsdB ; iedb = hi%IedB ; jsdb = hi%JsdB ; jedb = hi%JedB
512
513! This is where a control structure that is specific to this module would be allocated.
514 if (associated(cs)) then
515 call mom_error(warning, "register_restarts_dyn_unsplit_RK2 called with an associated "// &
516 "control structure.")
517 return
518 endif
519 allocate(cs)
520
521 alloc_(cs%diffu(isdb:iedb,jsd:jed,nz)) ; cs%diffu(:,:,:) = 0.0
522 alloc_(cs%diffv(isd:ied,jsdb:jedb,nz)) ; cs%diffv(:,:,:) = 0.0
523 alloc_(cs%CAu(isdb:iedb,jsd:jed,nz)) ; cs%CAu(:,:,:) = 0.0
524 alloc_(cs%CAv(isd:ied,jsdb:jedb,nz)) ; cs%CAv(:,:,:) = 0.0
525 alloc_(cs%PFu(isdb:iedb,jsd:jed,nz)) ; cs%PFu(:,:,:) = 0.0
526 alloc_(cs%PFv(isd:ied,jsdb:jedb,nz)) ; cs%PFv(:,:,:) = 0.0
527
528 thickness_units = get_thickness_units(gv)
529 flux_units = get_flux_units(gv)
530
531! No extra restart fields are needed with this time stepping scheme.
532
534
535!> Initialize parameters and allocate memory associated with the unsplit RK2 dynamics module.
536subroutine initialize_dyn_unsplit_rk2(u, v, h, tv, Time, G, GV, US, param_file, diag, CS, &
537 Accel_diag, Cont_diag, MIS, &
538 OBC, update_OBC_CSp, ALE_CSp, set_visc, &
539 visc, dirs, ntrunc, cont_stencil, dyn_h_stencil)
540 type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure.
541 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
542 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
543 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: u !< The zonal velocity [L T-1 ~> m s-1].
544 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(inout) :: v !< The meridional velocity [L T-1 ~> m s-1].
545 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Layer thicknesses [H ~> m or kg m-2]
546 type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic type
547 type(time_type), target, intent(in) :: time !< The current model time.
548 type(param_file_type), intent(in) :: param_file !< A structure to parse
549 !! for run-time parameters.
550 type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to
551 !! regulate diagnostic output.
552 type(mom_dyn_unsplit_rk2_cs), pointer :: cs !< The control structure set up
553 !! by initialize_dyn_unsplit_RK2.
554 type(accel_diag_ptrs), target, intent(inout) :: accel_diag !< A set of pointers to the
555 !! various accelerations in the momentum equations, which can
556 !! be used for later derived diagnostics, like energy budgets.
557 type(cont_diag_ptrs), target, intent(inout) :: cont_diag !< A structure with pointers
558 !! to various terms in the
559 !! continuity equations.
560 type(ocean_internal_state), intent(inout) :: mis !< The "MOM6 Internal State"
561 !! structure, used to pass around pointers
562 !! to various arrays for diagnostic purposes.
563 type(ocean_obc_type), pointer :: obc !< If open boundary conditions
564 !! are used, this points to the ocean_OBC_type
565 !! that was set up in MOM_initialization.
566 type(update_obc_cs), pointer :: update_obc_csp !< If open boundary
567 !! condition updates are used, this points
568 !! to the appropriate control structure.
569 type(ale_cs), pointer :: ale_csp !< This points to the ALE
570 !! control structure.
571 type(set_visc_cs), target, intent(in) :: set_visc !< set visc control structure
572 type(vertvisc_type), intent(inout) :: visc !< A structure containing
573 !! vertical viscosities, bottom drag
574 !! viscosities, and related fields.
575 type(directories), intent(in) :: dirs !< A structure containing several
576 !! relevant directory paths.
577 integer, target, intent(inout) :: ntrunc !< A target for the variable
578 !! that records the number of times the
579 !! velocity is truncated (this should be 0).
580 integer, intent(out) :: cont_stencil !< The stencil for
581 !! thickness from the continuity solver.
582 integer, intent(out) :: dyn_h_stencil !< The stencil for
583 !! thickness for the dynamics based on the
584 !! continuity solver and Coriolis scheme.
585
586 ! This subroutine initializes all of the variables that are used by this
587 ! dynamic core, including diagnostics and the cpu clocks.
588
589 ! Local variables
590 character(len=40) :: mdl = "MOM_dynamics_unsplit_RK2" ! This module's name.
591 character(len=48) :: flux_units
592 ! This include declares and sets the variable "version".
593# include "version_variable.h"
594 integer :: isd, ied, jsd, jed, nz, isdb, iedb, jsdb, jedb
595 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = gv%ke
596 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
597
598 if (.not.associated(cs)) call mom_error(fatal, &
599 "initialize_dyn_unsplit_RK2 called with an unassociated control structure.")
600 if (cs%module_is_initialized) then
601 call mom_error(warning, "initialize_dyn_unsplit_RK2 called with a control "// &
602 "structure that has already been initialized.")
603 return
604 endif
605 cs%module_is_initialized = .true.
606
607 cs%diag => diag
608
609 call log_version(param_file, mdl, version, "")
610 call get_param(param_file, mdl, "BE", cs%be, &
611 "If SPLIT is true, BE determines the relative weighting "//&
612 "of a 2nd-order Runga-Kutta baroclinic time stepping "//&
613 "scheme (0.5) and a backward Euler scheme (1) that is "//&
614 "used for the Coriolis and inertial terms. BE may be "//&
615 "from 0.5 to 1, but instability may occur near 0.5. "//&
616 "BE is also applicable if SPLIT is false and USE_RK2 "//&
617 "is true.", units="nondim", default=0.6)
618 call get_param(param_file, mdl, "BEGW", cs%begw, &
619 "If SPLIT is true, BEGW is a number from 0 to 1 that "//&
620 "controls the extent to which the treatment of gravity "//&
621 "waves is forward-backward (0) or simulated backward "//&
622 "Euler (1). 0 is almost always used. "//&
623 "If SPLIT is false and USE_RK2 is true, BEGW can be "//&
624 "between 0 and 0.5 to damp gravity waves.", &
625 units="nondim", default=0.0)
626
627 call get_param(param_file, mdl, "UNSPLIT_DT_VISC_BUG", cs%dt_visc_bug, &
628 "If false, use the correct timestep in the viscous terms applied in the first "//&
629 "predictor step with the unsplit time stepping scheme, and in the calculation "//&
630 "of the turbulent mixed layer properties for viscosity with unsplit or "//&
631 "unsplit_RK2. If true, an older incorrect value is used.", &
632 default=.false.)
633
634 call get_param(param_file, mdl, "DEBUG", cs%debug, &
635 "If true, write out verbose debugging data.", &
636 default=.false., debuggingparam=.true.)
637 call get_param(param_file, mdl, "TIDES", cs%use_tides, &
638 "If true, apply tidal momentum forcing.", default=.false.)
639 call get_param(param_file, mdl, "CALCULATE_SAL", cs%calculate_SAL, &
640 "If true, calculate self-attraction and loading.", default=cs%use_tides)
641
642 allocate(cs%taux_bot(isdb:iedb,jsd:jed), source=0.0)
643 allocate(cs%tauy_bot(isd:ied,jsdb:jedb), source=0.0)
644
645 mis%diffu => cs%diffu ; mis%diffv => cs%diffv
646 mis%PFu => cs%PFu ; mis%PFv => cs%PFv
647 mis%CAu => cs%CAu ; mis%CAv => cs%CAv
648
649 cs%ADp => accel_diag ; cs%CDp => cont_diag
650 accel_diag%diffu => cs%diffu ; accel_diag%diffv => cs%diffv
651 accel_diag%PFu => cs%PFu ; accel_diag%PFv => cs%PFv
652 accel_diag%CAu => cs%CAu ; accel_diag%CAv => cs%CAv
653
654 call continuity_init(time, g, gv, us, param_file, diag, cs%continuity_CSp, cs%OBC)
655 cont_stencil = continuity_stencil(cs%continuity_CSp)
656 call coriolisadv_init(time, g, gv, us, param_file, diag, cs%ADp, cs%CoriolisAdv)
657 dyn_h_stencil = max(cont_stencil, coriolisadv_stencil(cs%CoriolisAdv))
658 if (cs%calculate_SAL) call sal_init(h, tv, g, gv, us, param_file, cs%SAL_CSp)
659 if (cs%use_tides) call tidal_forcing_init(time, g, us, param_file, cs%tides_CSp)
660 call pressureforce_init(time, g, gv, us, param_file, diag, cs%PressureForce_CSp, cs%ADp, &
661 cs%SAL_CSp, cs%tides_CSp)
662 call hor_visc_init(time, g, gv, us, param_file, diag, cs%hor_visc)
663 call vertvisc_init(mis, time, g, gv, us, param_file, diag, cs%ADp, dirs, &
664 ntrunc, cs%vertvisc_CSp)
665 cs%set_visc_CSp => set_visc
666
667 if (associated(ale_csp)) cs%ALE_CSp => ale_csp
668 if (associated(obc)) cs%OBC => obc
669
670 flux_units = get_flux_units(gv)
671 cs%id_uh = register_diag_field('ocean_model', 'uh', diag%axesCuL, time, &
672 'Zonal Thickness Flux', flux_units, conversion=gv%H_to_MKS*us%L_to_m**2*us%s_to_T, &
673 y_cell_method='sum', v_extensive=.true.)
674 cs%id_vh = register_diag_field('ocean_model', 'vh', diag%axesCvL, time, &
675 'Meridional Thickness Flux', flux_units, conversion=gv%H_to_MKS*us%L_to_m**2*us%s_to_T, &
676 x_cell_method='sum', v_extensive=.true.)
677 cs%id_CAu = register_diag_field('ocean_model', 'CAu', diag%axesCuL, time, &
678 'Zonal Coriolis and Advective Acceleration', 'meter second-2', conversion=us%L_T2_to_m_s2)
679 cs%id_CAv = register_diag_field('ocean_model', 'CAv', diag%axesCvL, time, &
680 'Meridional Coriolis and Advective Acceleration', 'meter second-2', conversion=us%L_T2_to_m_s2)
681 cs%id_PFu = register_diag_field('ocean_model', 'PFu', diag%axesCuL, time, &
682 'Zonal Pressure Force Acceleration', 'meter second-2', conversion=us%L_T2_to_m_s2)
683 cs%id_PFv = register_diag_field('ocean_model', 'PFv', diag%axesCvL, time, &
684 'Meridional Pressure Force Acceleration', 'meter second-2', conversion=us%L_T2_to_m_s2)
685 cs%id_ueffA = register_diag_field('ocean_model', 'ueffA', diag%axesCuL, time, &
686 'Effective U-Face Area', 'm^2', conversion=gv%H_to_m*us%L_to_m, &
687 y_cell_method='sum', v_extensive=.true.)
688 cs%id_veffA = register_diag_field('ocean_model', 'veffA', diag%axesCvL, time, &
689 'Effective V-Face Area', 'm^2', conversion=gv%H_to_m*us%L_to_m, &
690 x_cell_method='sum', v_extensive=.true.)
691
692 id_clock_cor = cpu_clock_id('(Ocean Coriolis & mom advection)', grain=clock_module)
693 id_clock_continuity = cpu_clock_id('(Ocean continuity equation)', grain=clock_module)
694 id_clock_pres = cpu_clock_id('(Ocean pressure force)', grain=clock_module)
695 id_clock_vertvisc = cpu_clock_id('(Ocean vertical viscosity)', grain=clock_module)
696 id_clock_horvisc = cpu_clock_id('(Ocean horizontal viscosity)', grain=clock_module)
697 id_clock_mom_update = cpu_clock_id('(Ocean momentum increments)', grain=clock_module)
698 id_clock_pass = cpu_clock_id('(Ocean message passing)', grain=clock_module)
699 id_clock_pass_init = cpu_clock_id('(Ocean init message passing)', grain=clock_routine)
700
701end subroutine initialize_dyn_unsplit_rk2
702
703!> Clean up and deallocate memory associated with the dyn_unsplit_RK2 module.
704subroutine end_dyn_unsplit_rk2(CS)
705 type(mom_dyn_unsplit_rk2_cs), pointer :: cs !< dyn_unsplit_RK2 control structure that
706 !! will be deallocated in this subroutine.
707
708 dealloc_(cs%diffu) ; dealloc_(cs%diffv)
709 dealloc_(cs%CAu) ; dealloc_(cs%CAv)
710 dealloc_(cs%PFu) ; dealloc_(cs%PFv)
711
712 if (cs%calculate_SAL) call sal_end(cs%SAL_CSp)
713 if (cs%use_tides) call tidal_forcing_end(cs%tides_CSp)
714
715 deallocate(cs)
716end subroutine end_dyn_unsplit_rk2
717