Ice edge velocities for rheology and advection

[JFLemieux73]

This issue follows from a discussion with Bill Lipscomb.

a) Undefined (equal to 0) velocities are currently used at the ice edge for the calculation of strain rates. For this C-grid this affects the rheology term and diagnostics (div + shear). For the B-grid, this affects diagnostics (div + shear)...to be verified.

b) Does it also affect the U point u,v (for remap) obtained from interpolating uvelE and vvelN? To be verified.

There are two possible approaches:

1. Use masking logic to avoid spurious large gradients at cell edges where velocity components are not defined. 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.

2. Extend the mask where u and v are calculated. For these additional u and v, P=0 such that the freedrift velocities would be obtained. The momentum balance would be aice tauair = aice tauwater. A minimum value for aice might be required to avoid a division by zero in stepu (or step_vel).

Maybe the second approach is better as it deals with a) and b).

[dabail10]

So, we compute the B grid velocities from the C-grid velocities at the end of the subcycling:

call grid_average_X2Y('S', uvelE, 'E', uvel, 'U')
call grid_average_X2Y('S', vvelN, 'N', vvel, 'U')

These account for land points, but I'm positive that a valid ocean point with no ice would be averaged in. We do not use icetmask or iceumask in anyway here. I'm going to quick code something up and rerun my 5-year C-grid runs.

[dabail10]

Didn't really do much. The results are fairly similar.

[JFLemieux73]

I will get back to this soon. The calculation of uN and vE described in issue #976 could be useful here.

[JFLemieux73]

Here are a few free drift tests (P*=0).

kridge=-1
ktransport=-1.
ice_data_type = 'block'
ice_data_conc = 'c1'

We look at steady-state solution after a few days. We start with:

atm_data_type = 'uniform_east'

[JFLemieux73]

Here is uvel (B-grid):

This is good. vvel=0 and uvel is equal to the freedrift velocity everywhere there is ice.

[JFLemieux73]

The deformations (diagnostics) do not provide exactly the right information. Obviously the block is not deforming but here is the shear (B-grid again):

This is due to the zero velocities at the ice edge.

[JFLemieux73]

The divergence also exhibits non-zero values:

[JFLemieux73]

The C-grid solution also has uvelE=free drift value everywhere there is ice (not shown) and vvelN=0. I think that we don't see issues at the edges because vvelN=0. I will test later an oblique wind.

The divergence C-grid solution is like the B-grid one shown above (not shown). The shear is slightly different with kind of two grid cells with large values at the north and south edges (not sure why):

[JFLemieux73]

Now with oblique winds: atm_data_type = 'uniform_northeast'

B-grid: uvel=vvel everywhere there is ice (not shown). Similar issues with deformations. Here is the shear:

[JFLemieux73]

Now for the C-grid. Because of the problem discussed in issue #976 (#976 (comment)), the calculation of uN and vE close to the edges leads to non freedrift values at the edges. Here is uvelE:

Note that I changed the range and colormap to better see the different values at the edge.

[JFLemieux73]

Here is the shear for the C-grid: