Principal stress calculation does not work for C grid

[JFLemieux73]

The problem is that stress12T is used by principal_stress subroutine but it is never calculated (it is zero) for the C grid.

This affects sig1 and sig2 (as diagnotics). sigP is ok.

[apcraig]

I can help fix this. Can anyone provide guidance?

[eclare108213]

It looks like you could add a call in ice_history.F90 to get stress12T from wherever it's located on the C-grid. That might be the two edges, which would require an average? This is just a diagnostic.

[apcraig]

On the C grid, it looks like we compute stress12U. I assume that's on the U point. Should I just do an "average" of stress12U to get stress12T? Or is there a computation I can do directly on the T grid point?

[JFLemieux73]

The simplest and cheapest approach is to average the stress12U but I think the best approach (slightly more expensive...no issue for scaling) is to calculate stress12T inside the subcycling loop. You could just copy the calculation for stress12T in stressCD_T into stressC_T.

If we choose the averaging approach, it would be good to first compare the two approaches (i.e., plot the states of stress with ellipse).

[apcraig]

Thanks @JFLemieux73, I now see how to copy the stress12T calc into stressCD_T. I agree that's better. More soon.

[apcraig]

In stressC_T, we do

         shearTsqr = (shearU(i  ,j  )**2 * uarea(i  ,j  )  &
                    + shearU(i  ,j-1)**2 * uarea(i  ,j-1)  &
                    + shearU(i-1,j-1)**2 * uarea(i-1,j-1)  &
                    + shearU(i-1,j  )**2 * uarea(i-1,j  )) &
                    / (uarea(i,j)+uarea(i,j-1)+uarea(i-1,j-1)+uarea(i-1,j))

         DeltaT = sqrt(divT(i,j)**2 + e_factor*(tensionT(i,j)**2 + shearTsqr))

If we add the stress12T calculation, it needs shearT,

         stress12T(i,j) = (stress12T(i,j)*(c1-arlx1i*revp) &
                          + arlx1i*p5*etax2T(i,j)*shearT(i,j)) * denom1

should this be computed via

shearT = sqrt(shearTsqr)

which would always be positive (is that OK) or should we compute shearT via

      call strain_rates_T (nx_block   , ny_block     , &
                           icellT     ,                &
                           indxTi  (:), indxTj  (:)  , &
                           uvelE (:,:), vvelE   (:,:), &
                           uvelN (:,:), vvelN   (:,:), &
                           dxN   (:,:), dyE     (:,:), &
                           dxT   (:,:), dyT     (:,:), &
                           divT  (:,:), tensionT(:,:), &
                           shearT(:,:), DeltaT  (:,:))

the same way as in stressCD_T. At the present time, divT and tensionT are calculated in strain_rates_T in stressC_T. But we compute shearTsqr and DeltaT from shearTsqr locally as above.

If we compute shearT and deltaT via strain_rates_T in stressC_T, then stressC_T is identical to stressCD_T. Is that correct?

[JFLemieux73]

Hum...more complicated than I thought. I think we need the correct sign for shearT. So sqrt(shearT) does not work.

The problem if we want to calculate in stran_rates_T is that we don't have vvelE and uvelN to easily compute dvdx and dudy at the T point. Let me look a bit more closely.

[apcraig]

vvelE and uvelN are computed at the end of each subcycle. They are just a 4 pt average of vvelN and uvelE respectively. I don't know if that's good enough or not.

shearTsqr is really just an average of shearU**2. Maybe it's reasonable/correct to compute stress12T as an average of stress12U?

[JFLemieux73]

Or in stressC_T we also calculate shearT from the shearU values (similar to shearTsqr).

shearT=(shearU(i,j)* uarea(i ,j ) + ...

That's the method I would suggest. If we go with stress12T as an average of stress12U it would be good again to compare the states of stress with the other method.

[apcraig]

Just to clarify,

 shearT = (shearU(i  ,j  ) * uarea(i  ,j  )  &
                    + shearU(i  ,j-1) * uarea(i  ,j-1)  &
                    + shearU(i-1,j-1) * uarea(i-1,j-1)  &
                    + shearU(i-1,j  ) * uarea(i-1,j  )) &
                    / (uarea(i,j)+uarea(i,j-1)+uarea(i-1,j-1)+uarea(i-1,j))

in stressC_T? And shearTsqr /= shearT**2, that's OK, right?

[dabail10]

Since this is just a history diagnostic, I believe we should just average stress12U.

[JFLemieux73]

Yes in stressC_T.

And shearTsqr /= shearT**2, that's OK, right? ...you are right

Dave you are probably right but I would test both...I just want to make sure the average of the stress12U leads to stresses on and inside the yield curve.

[apcraig]

OK, I'll try both. I'll run C-grid cases with gx3 for 5 days. Is that enough? What do I need to output and plot?

[JFLemieux73]

@phil-blain you have some python code to plot the states of stress on the ellipse, right? Could you please make it available to Tony? Merci!

[apcraig]

I can make some elliptical plots, remind me what the two fields are I'm plotting. Is it sig1 vs sig2?

[phil-blain]

unfortunately I do not (Amélie had such a script which I did not copy in time before her account was deleted!)

[eclare108213]

I have an ancient one written in IDL - let me dredge it up

[eclare108213]

ellipse.pro.txt
This is ancient indeed. Are you familiar with IDL? Let me know if any of it is not understandable. Basically it defines and plots an ellipse with x and y axes, then reads in sig1 and sig2 and scatter-plots them on top.
The vectors xa and xb = [1, 2, ..., 300]/100 - 2.
imt, jmt are the array dimensions from the netcdf file (so 100, 116 for the gx3 grid).
The resulting plot will look something like this:

[apcraig]

OK, I have some plots. sig1 vs sig2.

• B grid
• C grid with shearT = avg(shearU)
• C grid with stress12T = avg(stress12U)

The shearT average which then goes into a stress12T computation is better. Its on/inside the ellipse. There are also points I don't show that are way off the domain of the plot for the avg(stress12U) case. So we should go with avg(shearU). It will be a little more expensive, but it's better and based on one test, looks like it won't increase cost by more than 0.5%.

One other note, the B-grid has 5x more points than the other 2 plots, so please ignore the density.

[JFLemieux73]

Thanks Tony. It's good to see that the points are on and inside the ellipse. Could you please try it with avg_strength instead of avg_zeta?

[apcraig]

One more plot

• C grid with shearT = avg(shearU) with visc_method=avg_strength

[JFLemieux73]

Great! Thanks again Tony. Having the stresses inside and on the ellipse is a good indication that what we coded is ok.