Generalize tracer_vertdiff and triDiagTs

[ashao]

tridiagTS is a slightly more optimized version of tracer_vertdiff (used for passive tracers)under the assumption that temperature and salinity have the same diapycnal diffusivity. This however leads to two similar pieces of code that must be maintained. Additionally, it has been suggested that the ea and eb input fields to tracer_vertdiff be replaced with Kd. For ALE cases, this change is simple since ea(k+1) = eb(k), however, this does not work for layered-configurations where in general ea(k+1) != eb(k) (due to diapycnal mass transport).

For full generality, I suggest the following:

1. Deprecate triDiagTS in favor of tracer_vertdiff, trading efficiency for generality
2. Make ea, eb, and Kd optional arguments to tracer_vertdiff and add two 3D arrays ea_loc, eb_loc
3. Perform the following logical tests at the top of tracer_vertdiff

• if ea and eb are both present: copy them to ea_loc and eb_loc
• if Kd is present: use it to calculate ea_loc and eb_loc
• if none are present or ea, eb, and Kd are all provided, throw a fatal error

Any other suggestions or alternatives?

[Hallberg-NOAA]

Thanks for these thoughts, @ashao.

I don't like the idea of adding additional 3-d arrays and copying into them, because that gets expensive. My preferred solution would be to exploit the fact that in ALE mode, ea(i,j,k) = 2*dt*Kd(i,j,K)/(h(i,j,k-1)+h(i,j,k)) is already essentially the second-order accurate time-integrated piston velocity that would be derived from Kd. In other words, once you know Kd, dt, and h, we can provide an interface that converts these three fields to ea. eb can be made into an optional argument to tracer_vertdiff when in ALE mode, for which eb(i,j,k)=ea(i,j,k-1), with the only loss being that the tracers would not then be added with any net mass flux from below.

Note that ea(i,j,1) incorporates net mass fluxes through the free surface, and that we do want the tracer fluxes to be consistent with net mass gain or loss through the free surface.