CICE-Consortium / CICE

Development repository for the CICE sea-ice model
Other
60 stars 132 forks source link

Instability at the coast in idealized C-grid eastblock simulation. #976

Closed JFLemieux73 closed 1 month ago

JFLemieux73 commented 2 months ago

While investigating issue #970 there was a crash after a few days with a negative area. I think this is different then all the cases we have seen before with the EFA method as these were related to the non-uniform grid and involved realistic velocities.

In this case the grid is uniform with dx=dy=16 km. The velocities are however unrealistic with some values 10 times larger than freedrift. This could be related to the problem seen by @daveh150 (issue #941).

The exp is:

./cice.setup --case INSTAB -g gbox80 -m ppp6 -e intel -p 1x1 -s boxforcee,boxclosed,buildincremental

I setup in ice_in the following:

ice_data_type = 'eastblock' ice_data_conc = 'c1' ice_data_dist = 'uniform'

Note that seabed_stress = .false.

With uniform_west (block drift away from the coast) there is no crash. I found the problem because I created a new idealized wind: NWW: uatm = -atm_val, vatm = -atm_val

JFLemieux73 commented 2 months ago

atm_data_type = 'uniform_NWW'.

JFLemieux73 commented 2 months ago

Ok I can really simplify the problem. Even with Pstar=0 it still crashes.

JFLemieux73 commented 2 months ago

Here is vvelN after 3 days. There are unrealistic values at the north and east edge, especially the NE corner.

E1_vvelN_2005-01-03-00000

JFLemieux73 commented 2 months ago

Keeping Pstar=0, I also set kridge = -1 ktransport = -1

Here is vvelN after 3 days

E2_vvelN_2005-01-03-00000

Notice the diff range of values. I just want to show that vvelN is different at the ice edge. It should not. It should be equal to the freedrift velocity. This is certainly related to our treatment of uvelE,vvelN at the ice edge.

JFLemieux73 commented 2 months ago

Ok the problem is clearly not related to the remapping. It also crashed with upwind.

JFLemieux73 commented 2 months ago

Ok I am able to correct the problem with ktransport=kridge=-1 for the diff velocities at the edge. At a point i,j (N) on the north edde, uN is (uEij+uEi-1j+0+0)/4. If I just use the ice points to calc uN, i.e. uN=(uEij+uEi-1j)/2 and I do the same thing for vE I get the following:

E3_vvelN_2005-01-03-00000

It is the freedrift velocity everywhere.

JFLemieux73 commented 2 months ago

Using this new code, it still crashes with ktransport=kridge=1. Looks like there is another problem...My bet is that it is related to the concentration used to calculate the wind stress and the water stress. With all the simplifications done, at steady state the balance is only between these two terms.

JFLemieux73 commented 2 months ago

At steady-state aicetauair = aicetauw. My guess is that aice on the LHS .ne. aice on RHS. In stepu_C and stepv_C I forced aice on both sides to be equal to 1. The problem disappears...

JFLemieux73 commented 1 month ago

Ok I think that aicetauair is only calculated at the beginning...i.e. aice = aice_init (constant in time). When solving for uvelE,vvelN at steady-state we would have:

tauw = aice_init x tauair / aice

When aice is small this can lead to large unrealistic velocities.

JFLemieux73 commented 1 month ago

@apcraig am I right that in idealized experiments such as the one described above aicetauair is only calculated once in uniform_data_atm and held fixed for the rest of the simulations?

JFLemieux73 commented 1 month ago

I just checked and the problem also affects the B-grid. The model also crashes.

apcraig commented 1 month ago

@apcraig am I right that in idealized experiments such as the one described above aicetauair is only calculated once in uniform_data_atm and held fixed for the rest of the simulations?

I don't know the answer for sure. But it looks like uniform_data_atm is called as part of the initialization (init_forcing_atmo) and during the run (get_forcing_atmo). When you added the "uniform_NWW" option, did you add it to both of those subroutines?

JFLemieux73 commented 1 month ago

I guess I created the bug with my new wind stress! Sorry about that.

Good to see that the idealized test cases are ok.

I will also modify it in get_forcing_atmo. Thanks!

JFLemieux73 commented 1 month ago

I am closing this issue as it was not a real instability. The only thing interesting here is the alternative calculation of uN and vE as explained above. This could be useful for issue #713.