Add Option to Specify Tracer Advection Time Step

[theresa-morrison]

Introduction

MOM6 separates dynamic and thermodynamic, or diabatic, processes and users have the option to specify a dynamic time step (DT) and a thermodynamic time step (DT_THERM). The dynamic time step is limited by the CLF, butDT_THERM is not and therefore is often longer than DT. Tracer advection is done using accumulated fluxes from one or more dynamic steps, such that the tracers are advected for the length of a thermodynamic step.

However, in higher resolution regional domains we have found that a tracer advective CFL may be providing a practical limit on tracer advection and therefore DT_THERM. We also know that model results are depend on the choice of DT_THERM (such as in Issue #773).

This PR introduces a new time step for tracer advection which can be less than or equal to the time step for diabatic processes. The tracer advection time step is DT_TRACER_ADVECT and can be an integer multiple of DT; DT_THERM should be an integer multiple of DT_TRACER_ADVECT. If DT_THERM is greater than the coupling time step, DT_TRACER_ADVECT can also be greater than the coupling time step. If DIABATIC_FIRST is true, DT_TRACER_ADVECT must equal DT_THERM, this limitation will be addressed in the near future. By default, DT_TRACER_ADVECT equals DT_THERM, which recovers old answers. The comments describing DT_THERM and DT_TRACER_ADVECT have been modified to reflect this change. Because the thermodynamics and tracer advection routines were already called separately there were no changes to those routines.

Results in NWA12

The sensitivity to DT_THERM can be clearly observed in the regional Northwest Atlantic 1/12th degree domain. These experiments all use shorter time steps in the first year, with DT=400s, and DT_THERM=800s. After the first year, DT=600s and DT_THERM is set to 1800s, or 3600s. Using a tracer advection time step shorter than 3600s improved the Gulf Stream position with DT_THERM=3600, but more detailed analysis is needed. Here the position of the Gulf Stream is shown for the different combinations of timesteps.

Experiment 1: DT=600, DT_THERM=1800

Experiment 2: DT=600, DT_THERM=3600

Experiment 3: DT=600, DT_TRACER_ADVECT=1800, DT_THERM=3600

Experiment 4: DT=600, DT_TRACER_ADVECT=1200, DT_THERM=3600

Testing

This change was tested in the 0.25 degree Baltic domain for reproducibility across restarts for a combination of time step options:

• DT = 300, DT_TRACER_ADVECT = 600, DT_THERM = 1200, DT_coupled = 2400
• DT = 300, DT_TRACER_ADVECT = 600, DT_coupled = 1200, DT_THERM = 2400
• DT = 300, DT_coupled = 600, DT_TRACER_ADVECT = 1200, DT_THERM = 2400

Note that the new option TRADV_SPANS_COUPLING has been added, and should be used where appropriate. Because grid generation and remapping are done in step_MOM_thermo at the end of every thermodynamic time step, DT_TRACER_ADVECT must be less than DT_THERM.

[awallcraft]

I don't like DT_TADVECT as a name, I would call this DT_ADVECT or perhaps DT_TRADVECT

The description of DR_THERM should say "Ideally DT_THERM should be an integer multiple of DT and of DT_TADVECT "

It seems that if the coupling time step is 3600, for example, then DT_TADVECT=10000 DT_THERM=20000, both SPANS_COUPLING, would be silently converted to DT_TADVECT=7200 and DT_THERM=18000 and the later isn't an integer multiple of DT_TADVECT.

[Hallberg-NOAA]

Overall I really like this contribution, and I think that the figures in the conversation demonstrate what a valuable addition this is. However, I concur with the suggestion that the new variable should be called DT_TR_ADVECT instead of DT_TADVECT.

Also, I think that it should be possible to use tracer timesteps that are not the same with each step, so that in the case of Alan's example where DT_THERMO is reset to 18000 s but DT_TR_ADVECT becomes 7200 s, the actual sequence of tracer advection timesteps would be 7200 s, 7200 s, 3600 s within a single DT_THERMO of 18000 s. I am confident that the iterative tracer advection scheme can handle this, and the logic setting the tracer advection and thermodynamics timesteps should not be too bad. I would, however, be satisfied if this option were held for another PR, especially if for now we added some error handling to deal with the cases where DT_THERMO is not an integer multiple of DT_TR_ADVECT

[theresa-morrison]

Thanks Alan and Bob for the suggestions. I am happy to rename the time step to make things clearer.

[theresa-morrison]

I agree some of the reccomendations for the timesteps in their descriptions could be modified and the possible options for timesteps could be handled better.

The description of DR_THERM should say "Ideally DT_THERM should be an integer multiple of DT and of DT_TADVECT "

I can make this change as well.

Also, I think that it should be possible to use tracer timesteps that are not the same with each step, so that in the case of Alan's example where DT_THERMO is reset to 18000 s but DT_TR_ADVECT becomes 7200 s, the actual sequence of tracer advection timesteps would be 7200 s, 7200 s, 3600 s within a single DT_THERMO of 18000 s.

I believe the current logic would support this example if thermo spans coupling is false. In that case, if DT=3600, DT_THERMO= 18000s, and DT_TR_ADVECT=7200 then ntastep would be 2. When thermo spans coupling is false, the check depends on the step number n relative to ntastep and always happens on the last step. This means the tracer advection timestep would effectively be: 7200s, 7200s, 3600s in a single tracer timestep.

If thermo spans coupling is true, the elapsed time since the last advection step is checked. But I will add an additional requirement to do advection if the t_dyn_rel_thermo==DT_THERM.

[theresa-morrison]

I have updated the variable names in the code and the description of this PR. DT_TADVECT is now DT_TRACER_ADVECT, dt_tadvect is now dt_tr_adv, and tadvect_spans_coupling is now tradv_spans_coupling.

I added line 949:
if (CS%t_dyn_rel_thermo + 0.5*dt > dt_therm) do_advection = .true.
which should force a shorter advection step to occur if the time is close enough to the end of the thermodynamic loop.

I tested these changes for reproducibility across restarts in the 0.25 degree Baltic domains with the time steps listed above and two new options:

• DT = 300, DT_TRACER_ADVECT = 900, DT_coupled = 1800, DT_THERM = 2400
• DT = 300, DT_coupled = 600, DT_TRACER_ADVECT = 900, DT_THERM = 2400
  to test some uneven steps. The tests ran successfully and passed the restart tests.

[Hallberg-NOAA]

Thank you for this valuable new capability to set the tracer advection and thermodynamic timesteps separately.

[Hallberg-NOAA]

In the coupled_AM2_LM3_SIS2/Intersperse_ice_1deg test case, this new code is failing with what the intel compiler reports as an integer division by zero, perhaps at line 968 of MOM.F90.

[yichengt900]

@theresa-morrison, I conducted several runtime experiments for this new feature in both physics-only and OBGC-coupled modes using the NWA domains. The summarized results are presented in the table below. In the physics-only configuration, the runtime improvement was modest, achieving only a 4% reduction. However, in the OBGC-coupled mode, the runtime showed a significant reduction of approximately 17-19% (dachieved by doubling DT_THERM while keeping DT_TRACER_ADVECT close to the CTRL run). The results maintained ok Gulf Stream positioning and consistency across other OBGC tracers (not shown here).

Avg runtime (mins)

Mode	CTRL (DT=600s, DT_TRACER_ADVECT=DT_THERM=1800s)	DT=600s, DT_TRACER_ADVECT=1200s, DT_THERM=3600s	speedup
MOM6	93	89	-4%
MOM6-OBGC	235	191	-19%

[theresa-morrison]

In the coupled_AM2_LM3_SIS2/Intersperse_ice_1deg test case, this new code is failing with what the intel compiler reports as an integer division by zero, perhaps at line 968 of MOM.F90.

I was able to recreate this error in an interspersed Baltic test case. The solution I found was to add an additional check that do_thermo was true before checking if do_diabatic is true. If the order of step_MOM calls is being set by the combined ice ocean driver, the model should only consider doing a diabatic update when doing thermodynamics. One of the variables (ntstep) used in the do_diabatic check is not set if do_thermo is false, leading to a division by zero error.

[Hallberg-NOAA]

It now appears that there are cases where do_diabatic is not being initialized, but it is being used on line 972. Logically the value of do_diabatic does not matter, but the use of an uninitialized variable will be flagged when certain debugging capabilities are enabled. Please consider adding else ; do_diabatic = .false. just before line 971.

There are changes being implemented and a failed test so I'm removing the approval to avoid confusion