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