New method for calculating C-grid viscosities

[JFLemieux73]

C-grid viscosities are currently calculated at the T point following Bouillon et al. 2013 (see also Kimmritz et al. 2016).
At the U point, eta can be calculated following Kimmritz et al. 2016 C1 approach (avg_zeta) and C2 approach (avg_strength).

Bill Lipscomb suggests a new method described below:

This is a suggested method for computing and averaging the various stress, strain-rate, and strength terms during the EVP subcycling. The goals are:
Use masking logic to avoid spurious large gradients at cell edges where velocity components are not defined.
Minimize the averaging back and forth between T and U points.
If possible, average only quantities that are grid-invariant (i.e., unchanged under rotation). Recall that DD is a grid invariant, as is sqrt(DT2 + DS2), and also Delta = sqrt[DD2 + (1/e2)(DT2 + DS2)]. Scalars like P are also invariant.

Here are the proposed steps:
Compute P at T points in the usual way. Set P = 0 in ice-free ocean cells.
Average P from T points to U points. Include ice-free cells (P = 0) in the average.
Compute du/dx, dv/dy, DD and DT at T points in the usual way.
Compute du/dy, dv/dx and DS at U points in the usual way, except:
If one of the two adjacent edges is masked out (because it sits between two ice-free ocean cells with indeterminate speed), set the gradient (either du/dy or dv/dx) to zero. If one of the edges is land-based with a speed of zero, it's OK to include the zero value in the average.
Compute Delta in each quadrant of each ice-filled cell, using DD and DT from the T point, and using DS from the U point. In a given cell, each quadrant will have a different Delta, since the four vertices each have a different DS. Note that Delta is undefined in ice-free cells.
At each T point and each U point, compute Delta as the average of the values in the 4 neighboring quadrants. When averaging, do not include ice-free cells where Delta is undefined.
Compute replacement pressure, zeta, and eta at T points based on P and Delta at T points.
Compute eta at U points based on P and Delta at U points.
Compute stressp at T points based on zetaT, divergT and rep_prsT.
Compute stressm at T points based on etaT and tensionT.
Compute stress12 at U points based on etaU and shearU.
Compute div_stress_Ex and div_stress_Ny as in the current code, and so on till the end.