From 4bf730a7a8ec9fb2fedb00476cdd878adb944571 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Thu, 23 Sep 2021 19:02:28 -0400 Subject: [PATCH] +(*)Add ALTERNATE_FIRST_DIRECTION Added a new runtime parameter, ALTERNATE_FIRST_DIRECTION, to cause the first direction in the directionally split operators (such as the continuity solver, the barotropic solver, and the tracer advection) to alternate with every dynamic timestep. This includes changes to record this direction in the restart file, so that the model will reproduce across restarts. In addition a previously uninitialized logical in offline_advection_layer is now being initialized similarly to other similar variables in the MOM_offline_main module. By default, all answers in the MOM6-examples are bitwise identical, but there is a new runtime parameter and a new non-mandatory field in the MOM6 restart files, and there might be changes in offline tracer advection runs. --- src/core/MOM.F90 | 39 +++++++++++++++++++++++++++++++-- src/tracer/MOM_offline_main.F90 | 1 + 2 files changed, 38 insertions(+), 2 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index eea888cd70..aa593cba7d 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -212,6 +212,8 @@ module MOM real :: t_dyn_rel_adv !< The time of the dynamics relative to tracer advection and lateral mixing !! [T ~> s], or equivalently the elapsed time since advectively updating the !! tracers. t_dyn_rel_adv is invariably positive and may span multiple coupling timesteps. + integer :: n_dyn_steps_in_adv !< The number of dynamics time steps that contributed to uhtr + !! and vhtr since the last time tracer advection occured. real :: t_dyn_rel_thermo !< The time of the dynamics relative to diabatic processes and remapping !! [T ~> s]. t_dyn_rel_thermo can be negative or positive depending on whether !! the diabatic processes are applied before or after the dynamics and may span @@ -233,6 +235,9 @@ module MOM logical :: use_ALE_algorithm !< If true, use the ALE algorithm rather than layered !! isopycnal/stacked shallow water mode. This logical is set by calling the !! function useRegridding() from the MOM_regridding module. + logical :: alternate_first_direction !< If true, alternate whether the x- or y-direction + !! updates occur first in directionally split parts of the calculation. + real :: first_dir_restart = -1.0 !< A real copy of G%first_direction for use in restart files logical :: offline_tracer_mode = .false. !< If true, step_offline() is called instead of step_MOM(). !! This is intended for running MOM6 in offline tracer mode @@ -1137,6 +1142,11 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, & ! Advance the dynamics time by dt. CS%t_dyn_rel_adv = CS%t_dyn_rel_adv + dt + CS%n_dyn_steps_in_adv = CS%n_dyn_steps_in_adv + 1 + if (CS%alternate_first_direction) then + call set_first_direction(G, MODULO(G%first_direction+1,2)) + CS%first_dir_restart = real(G%first_direction) + endif CS%t_dyn_rel_thermo = CS%t_dyn_rel_thermo + dt if (abs(CS%t_dyn_rel_thermo) < 1e-6*dt) CS%t_dyn_rel_thermo = 0.0 CS%t_dyn_rel_diag = CS%t_dyn_rel_diag + dt @@ -1168,6 +1178,7 @@ subroutine step_MOM_tracer_dyn(CS, G, GV, US, h, Time_local) !! of the time step. type(group_pass_type) :: pass_T_S integer :: halo_sz ! The size of a halo where data must be valid. + logical :: x_first ! If true, advect tracers first in the x-direction, then y. logical :: showCallTree showCallTree = callTree_showQuery() @@ -1189,8 +1200,16 @@ subroutine step_MOM_tracer_dyn(CS, G, GV, US, h, Time_local) call cpu_clock_begin(id_clock_thermo) ; call cpu_clock_begin(id_clock_tracer) call enable_averages(CS%t_dyn_rel_adv, Time_local, CS%diag) + if (CS%alternate_first_direction) then + ! This calculation of the value of G%first_direction from the start of the accumulation of + ! mass transports for use by the tracers is the equivalent to adding 2*n_dyn_steps before + ! subtracting n_dyn_steps so that the mod will be taken of a non-negative number. + x_first = (MODULO(G%first_direction+CS%n_dyn_steps_in_adv,2) == 0) + else + x_first = (MODULO(G%first_direction,2) == 0) + endif call advect_tracer(h, CS%uhtr, CS%vhtr, CS%OBC, CS%t_dyn_rel_adv, G, GV, US, & - CS%tracer_adv_CSp, CS%tracer_Reg) + CS%tracer_adv_CSp, CS%tracer_Reg, x_first_in=x_first) call tracer_hordiff(h, CS%t_dyn_rel_adv, CS%MEKE, CS%VarMix, G, GV, US, & CS%tracer_diff_CSp, CS%tracer_Reg, CS%tv) if (showCallTree) call callTree_waypoint("finished tracer advection/diffusion (step_MOM)") @@ -1213,6 +1232,7 @@ subroutine step_MOM_tracer_dyn(CS, G, GV, US, h, Time_local) call cpu_clock_begin(id_clock_thermo) ; call cpu_clock_begin(id_clock_tracer) CS%uhtr(:,:,:) = 0.0 CS%vhtr(:,:,:) = 0.0 + CS%n_dyn_steps_in_adv = 0 CS%t_dyn_rel_adv = 0.0 call cpu_clock_end(id_clock_tracer) ; call cpu_clock_end(id_clock_thermo) @@ -1994,6 +2014,11 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & "in parts of the code that use directionally split "//& "updates, with even numbers (or 0) used for x- first "//& "and odd numbers used for y-first.", default=0) + call get_param(param_file, "MOM", "ALTERNATE_FIRST_DIRECTION", CS%alternate_first_direction, & + "If true, after every dynamic timestep alternate whether the x- or y- "//& + "direction updates occur first in directionally split parts of the calculation. "//& + "If this is true, FIRST_DIRECTION applies at the start of a new run or if "//& + "the next first direction can not be found in the restart file.", default=.false.) call get_param(param_file, "MOM", "CHECK_BAD_SURFACE_VALS", CS%check_bad_sfc_vals, & "If true, check the surface state for ridiculous values.", & @@ -2256,6 +2281,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & ALLOC_(CS%uhtr(IsdB:IedB,jsd:jed,nz)) ; CS%uhtr(:,:,:) = 0.0 ALLOC_(CS%vhtr(isd:ied,JsdB:JedB,nz)) ; CS%vhtr(:,:,:) = 0.0 CS%t_dyn_rel_adv = 0.0 ; CS%t_dyn_rel_thermo = 0.0 ; CS%t_dyn_rel_diag = 0.0 + CS%n_dyn_steps_in_adv = 0 if (debug_truncations) then allocate(CS%u_prev(IsdB:IedB,jsd:jed,nz)) ; CS%u_prev(:,:,:) = 0.0 @@ -2395,6 +2421,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & ! Set a few remaining fields that are specific to the ocean grid type. call set_first_direction(G, first_direction) + CS%first_dir_restart = real(G%first_direction) ! Allocate the auxiliary non-symmetric domain for debugging or I/O purposes. if (CS%debug .or. G%symmetric) then call clone_MOM_domain(G%Domain, G%Domain_aux, symmetric=.false.) @@ -2450,6 +2477,12 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & CS%tv%S => CS%S endif + ! Reset the first direction if it was found in a restart file. + if (CS%first_dir_restart > -0.5) & + call set_first_direction(G, NINT(CS%first_dir_restart)) + ! Store the first direction for the next time a restart file is written. + CS%first_dir_restart = real(G%first_direction) + call rotate_initial_state(u_in, v_in, h_in, T_in, S_in, use_temperature, & turns, CS%u, CS%v, CS%h, CS%T, CS%S) @@ -2720,7 +2753,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & call insert_offline_main( CS=CS%offline_CSp, ALE_CSp=CS%ALE_CSp, diabatic_CSp=CS%diabatic_CSp, & diag=CS%diag, OBC=CS%OBC, tracer_adv_CSp=CS%tracer_adv_CSp, & tracer_flow_CSp=CS%tracer_flow_CSp, tracer_Reg=CS%tracer_Reg, & - tv=CS%tv, x_before_y = (MOD(first_direction,2)==0), debug=CS%debug ) + tv=CS%tv, x_before_y=(MODULO(first_direction,2)==0), debug=CS%debug ) call register_diags_offline_transport(Time, CS%diag, CS%offline_CSp) endif @@ -3007,6 +3040,8 @@ subroutine set_restart_fields(GV, US, param_file, CS, restart_CSp) "Density unit conversion factor", "R m3 kg-1") call register_restart_field(US%J_kg_to_Q_restart, "J_kg_to_Q", .false., restart_CSp, & "Heat content unit conversion factor.", units="Q kg J-1") + call register_restart_field(CS%first_dir_restart, "First_direction", .false., restart_CSp, & + "Indicator of the first direction in split calculations.", "nondim") end subroutine set_restart_fields diff --git a/src/tracer/MOM_offline_main.F90 b/src/tracer/MOM_offline_main.F90 index 408120b4e5..5f1d59337e 100644 --- a/src/tracer/MOM_offline_main.F90 +++ b/src/tracer/MOM_offline_main.F90 @@ -895,6 +895,7 @@ subroutine offline_advection_layer(fluxes, Time_start, time_interval, CS, h_pre, IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB dt_iter = CS%US%s_to_T * time_interval / real(max(1, CS%num_off_iter)) + x_before_y = CS%x_before_y do iter=1,CS%num_off_iter