Refactor rheology (stress) computations

[phil-blain]

I went through the code and there is a lot of repetitions (copy-paste) between different computations regarding rheology. Here is a list of the computations by subroutines:

calc_zeta_dPr (ice_dyn_vp)

• computes zeta
• computes \partial_x P_R, \partial_y Pr (derivatives of the replacement pressure)
• computes stressp_[1-4] (only P_R part)
• then computes [cs]sigp[nsew12]
• stress_Pr -> used in calc_bvec

stress_vp (ice_dyn_vp)

• computes stresses from zeta
• stress[p,m,12]_[1-4]
• here stressp is the whole term

matvec (ice_dyn_vp)

• computes (partial stresses) from zeta
• stress[p,m,12]_[1-4]
• stressp is only "linear" part (linear in \partial^n_[x,y] u)
• computes [cs]sigp[nsew12]
• computes str

formDiag_step1 (ice_dyn_vp)

• computes strain_rates (with some coefficients = 0 or 1)
• computes stresses (stressp -> only "linear" part)
• computes Drheo (same computation as matvec, but special case for 8 corners...)

formDiag_step2 (ice_dyn_vp)

• similar to the second step of matvec, bu without ccb (i.e. only ccaimp)

stress (ice_dyn_evp)

• computes stress[p,m,12]_[1-4] -> different computations (EVP)
• computes [cs]sig[p,m,12][nsew] -> SAME as matvec
• computes str[pm]_tmp -> SAME as matvec

stress_eap (ice_dyn_eap)

• computes strain_rates
• computes deformations
• computes [cs]sig[p,m,12][nsew] -> SAME as matvec
• computes str[pm]_tmp -> SAME as matvec

We will refactor the code to unify these computations, but in a second step (after the Picard PR).