Modify uN and vE for for the computation of seabed stress coeff on the C-grid.

[JFLemieux73]

uN and vE for the computation of the Cb coeff are for the moment weighted averages of uE and vN. This is inconsistent with the C-grid approach in Lemieux et al. 2015. It should be:

|uN|=min(|uEij|, |uEij+1|, |uEi-1j|, |uEi-1j+1|), same idea for vE.

We should also create idealized test cases for grounding and see the impact of this change.

[JFLemieux73]

I modified stepu_C and stepv_C to calc |uN| and |vE| has explained above. This is in my branch GroundingModif.

I coded a simple grounding test with a high shoal (hwater = 5 m) at iglob= 74 (gbox80). The ice is initialize with:

ice_data_type = 'eastblock'
atm_data_type = 'uniform_east'
ice_data_conc = 'c1'

Here is hi after 2 weeks with the current code:

[JFLemieux73]

And here is hi with the new code:

[JFLemieux73]

The solutions are very similar. The larger block is drifting while the smaller rectangle is landfast as there is grounding at i=74.

[JFLemieux73]

Hum...I just realized that the new code does not give the same answer when changing the decomposition. To fix the problem the halo update

call dyn_haloUpdate (halo_info, halo_info_mask, &
field_loc_Eface, field_type_vector, &
uvelE)

should be moved between call stepu_C and call stepv_C.

[JFLemieux73]

The current code has:

do iblk = 1, nblocks
     call stepu_C..
     call stepv_C..
enddo
call dyn_haloUpdate (uvelE)
call dyn_haloUpdate (vvelN)

while the new code would have:

do iblk = 1, nblocks
    call stepu_C..
enddo
call dyn_haloUpdate (uvelE)
do iblk = 1, nblocks
     call stepv_C..
enddo
call dyn_haloUpdate (vvelN)

It is a bit more complicated than anticipated...

[JFLemieux73]

I just tested the modification above and the answers are BFB when changing the decomposition.

[JFLemieux73]

I am taking a step back and will look at a simpler problem. There is shoal at i=74 and a sea ice cover from i=74 up to the coast on the right (east). I call it minieastblock. The wind is uniform_west. I start with the current code. Here is hi after 3 days.

Note that the range is between 0.9 and 1m...Results look good; it is landfast, velocities are small (not shown). There is a little bit of 'diffusion' at the top and bottom but overall the current code does a good job.

[JFLemieux73]

And here is with the new code:

Very similar...

[JFLemieux73]

Now with oblique winds from the north-east (current code).

It is landfast except at the open (i.e. no shoal) south edge. The ice probably fails in tension-shear at this location. Note that the range is now 0-1m.

[JFLemieux73]

Now with the new code:

Again results are quite similar. Why should we modify the code then?

[JFLemieux73]

I think this was important with the McGill model because the ice edge condition (freedrift) is different than the current one in CICE (for the C-grid), that is uvelE=vvelN=0 when there is no ice (or aice<amin).

[JFLemieux73]

Hi @eclare108213, I am thinking of leaving the code as is. The current code does a good job for the seabed stress (at least in these idealized experiments). The calculation of uN and vE for the seabed stress is simpler and than in Lemieux et al. 2015. I guess the approach of Lemieux et al. 2015 is not required because the treatment at the ice edge is different (see last comment above). Do you have any concerns?

[eclare108213]

@JFLemieux73, two thoughts. First, you pointed out in the C-grid paper that the implementation is different from your 2015 paper and would be fixed. I'm less concerned about that than the possibility that someone implements free drift for the ice edge and then has a problem with this bit of code. The simulations are very similar in these tests, but is the McGill version better in a physical sense? I've always thought that we ought to be using free drift for the ice edge, but (correct me if I'm wrong) so far no one has made a convincing argument for doing so. So I lean toward changing the code if it's more physical and will head off future problems, but I'm also okay with leaving the code as-is.

[JFLemieux73]

If you look at this part of issue #976

#976 (comment)

you can see that the abrupt ice edge leads to unrealistic velocity values on the contour. I found a fix with uN=( uEij + uEi-1j )/2 but the freedrift velocity would also work with uN=( uEij + uEi-1j + ufreedrift + ufreedrift )/4