Coupled MOM6/SIS2: problem in ATM -> ICE/OCN coupling

[npoumaere]

Hi everyone,

Hope you're all doing well.

I am trying to set up a coupled MOM6/SIS2 Southern Ocean channel model.
I started from the Baltic example in ice_ocean_SIS2 folder, and modified the input files to match my case.

The configuration is as follows:
- Extent: -70S -30S latitude, 0E 50E longitude, zonally reentrant
- Coarse 1 degree resolution (I want to get the coupling right before stepping up the resolution)
- Sponge at northern boundary on Theta and S. Free slip at southern boundary.
- Monthly varying surface forcings of wind (direct wind stress or 10m wind velocity, see problem below), air surface temperature and specific humidity, downward shortwave and longwave radiation, rain and snow precipitation.

Forcings are time varying, zonally invariant profiles broadly derived from SOSE and ERA5 datasets.
The model runs for at least 10 years without crashing.

PROBLEM:

There seems to be something wrong happening in the coupling, because:
- The ocean -> ice heat flux is of the order of +1000/+5000 W/m2, and the SSTs get colder and colder regardless of the imposed surface air temperature.
- Despite this flux, sea ice grows unbounded both in thickness (reaching 5m at southern boundary over 10 years) and northern extent (reaching -35degS in 10 years). There is a seasonal melting/freezing cycle though.
- Currents are coherent with mechanical forcing when direct wind stress is applied to the ICE component in data_table, but switching to imposed 10m wind speed in ATM component leads to no wind stress at all and thus no circulation.

So my guess is that there is no coupling between the dummy atmosphere component and the OCN/ICE components, hence no sensible heat flux and wind stress.
The tremendous OCN->ICE heat flux hints at another problem, but I can't see what.

I tried several things:
- changing the specification of the *tile1X* files by imposing some land cells, in case that would be needed, but no change in the results. For now I follow advice from this thread to generate *tile1X* with ncells=0 for ATM/LND and LND/OCN files.
- specifying direct heat flux on the ICE component (t_flux in data_table), but that leads to crash after one year.
- Switching on and off the "do_flux" option in coupler_nml (it is off in Baltic, on in my config), but switching it off leads to a crash.

It feels like I'm missing something obvious, but I don't see what that could be.
If someone has some insight into that kind of problem I'd be glad to hear any clue or possible avenue to solve that.

Attached is a compressed archive containing all input files, including surface forcing. (the unpacked folder weighs 26M, you can run "tar -xzvf SO_FILES.tar.gz" to unpack it)

All the best,
Nelson

 

[npoumaere]

Hi everyone,

I think I found where the problem came from.

Basically, I had set the gustiness field to zero everywhere, and my 10m wind speed field has some regions of zero values. My guess is that having both the wind speed and gustiness being zero leads to an error somewhere in the code (can't say exactly where yet) which ended up in crashes, so I had to resort to directly specifying the wind stress to the ICE component (the u_flux and v_flux in the data_table). But the computation of the sensible heat flux (related to the specified atmosphere bottom temperature) seems to be directly impacted by the wind speed field (no wind ~= no sensible heat flux), hence discrepancies and uncontrolled sea ice growth.
Setting the gustiness to a small constant value (0.001) in data_table seems to be doing the trick: allowing to specify 10m wind speed and with coherent wind stress and sensible heat flux being computed from it.
Setting gustiness to zero (in data_table) leads to the following error message: "NaN detected: u Before steps forces%tau[xy]".

The first results from a 10-year run are encouraging, with expected sea ice growth/melt cycle and no excessive cooling and freshening of the basin surface.

Still, one strange thing is that the specification of the *tile1X* files has no impact on the computation at all: for the ATM/LND and LND/OCN files, setting ncells=0 or ncells=NIGLOBAL*NJGLOBAL with xgrid_area = Ah[i,j] doesn't change the results in the slightest (same fields computed down to machine precision), which is a bit confusing.
I wonder if someone knows about this kind of behaviour?
(Note that when setting ncells=0, I have to subsequentially run " ncks --fix_rec_dmn ncells " to make the *tile1X* files compatible with the code, which actually sets ncells to 1.)

Attached is the updated config.

All the best,
Nelson