Here is a question about the computation of W before
time-stepping tracer fields. Let's consider a b/c time step, from time
n to n+1. To time-step tracers, mass fluxes are needed at the
intermediate time n+1/2 (is my interpretation correct?).
Horizontal fluxes are updated in step3d_uv, and vertical fluxes in omega,
(when omega is called for the 2nd time in main3d).
Does it look like that during this second call the horizontal fluxes
(Huom, Hvon) correspond to time n+1/2 and layer heights (z_w)
to a different time, n+1?
Is that true? If so, is it ok to use Huon and z_w at different times?
Below are some details that may provide some background info.
As far as I understand the code, after the b/t loop is done, we obtain:
(1) the filtered SSH at time n+1 (Zt_avg1);
(2) the filtered depth-integrated volume flux at time n+1
(DU_avg1, DV_avg1);
(3) the filtered depth-integrated volume flux at time n+1/2
(DU_avg2, DV_avg2).
(4) updated s-layer locations z_w, z_r and layer depths Hz
at time n+1 (those are obtained using set_depth and Zt_avg1).
Next, step3d_uv updates u and v (at time n+1) and adjusts their vert.
integral to be consistent with DU_avg1, DV_avg1, and Hz.
This subroutine (step3d_uv) also computes mass fluxes Huon, Hvom at
time n+1/2, and adjusts their vertical integral using DU_avg2, DV_avg2.
These are used, together w/ z_w [see (4)], to compute W (omega.F).
Then, Huon, Hvom, and W are used to time-step tracers.