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

 

[ruyin Wang]

Hi Nelson,

How do you change the data_table file?

Thanks for answering

 

[npoumaere]

Hi Ruyin,

Could you please be a bit more specific? My data_table is available in the above attached SO_FILES_UPDATE.tar.gz if that can help.

All the best,
Nelson

 

[ruyin Wang]

Hi Ruyin,

Could you please be a bit more specific? My data_table is available in the above attached SO_FILES_UPDATE.tar.gz if that can help.

All the best,
Nelson

Thank you for your response. I have noticed a couple of additional questions:

1. It seems that you have replaced the forcing files in the data_table with your own datasets. Is that correct? The original configuration appears to use the CORE NYF forcing.
2. A slightly unrelated question: I noticed that you modified NIGLOBAL = 50 and NJGLOBAL = 68 in the MOM_input. In one of my experiments, I also tried modifying these parameters using MOM_override, but it seems to have caused some grid inconsistencies. May I ask whether you also modified files such as ocean_hgrid.nc, topog.nc, etc., accordingly?

Thank you very much for your help ！！

 

[npoumaere]

1. Yes, you can tell MOM6 to use the forcing dataset you want, whether you generated it yourself or not, through the data_table. The dimensions of the forcing dataset simply need to match the ones you specify in MOM_input/override

2. Yes, the NIGLOBAL and NJGLOBAL parameters need to match the zonal and meridional dimensions in ocean_hgrid.nc, topog.nc, etc. If you want to change the resolution i.e. change NIGLOBAL/NJGLOBAL, you need to modify ocean_hgrid.nc accordingly, along with all other geometry/forcing files.

Note that ocean_hgrid.nc is a "mosaic" file, and thus specifies the "supergrid", that is both the tracer (or thickness) points and the velocity points (see Discrete Horizontal and Vertical Grids — MOM6 0.2a3 documentation). So for example the number of points in the zonal direction in this file needs to be 2*NIGLOBAL + 1.

Regards,
Nelson

 

[ruyin Wang]

1. Yes, you can tell MOM6 to use the forcing dataset you want, whether you generated it yourself or not, through the data_table. The dimensions of the forcing dataset simply need to match the ones you specify in MOM_input/override

2. Yes, the NIGLOBAL and NJGLOBAL parameters need to match the zonal and meridional dimensions in ocean_hgrid.nc, topog.nc, etc. If you want to change the resolution i.e. change NIGLOBAL/NJGLOBAL, you need to modify ocean_hgrid.nc accordingly, along with all other geometry/forcing files.

Note that ocean_hgrid.nc is a "mosaic" file, and thus specifies the "supergrid", that is both the tracer (or thickness) points and the velocity points (see Discrete Horizontal and Vertical Grids — MOM6 0.2a3 documentation). So for example the number of points in the zonal direction in this file needs to be 2*NIGLOBAL + 1.

Regards,
Nelson

Thanks so much for your help, and I have one more question that I see in your input.nml file, you set the "calendar = 'thirty_day'". What is this calender mean and what's different between the "NOLEAP"

Thanks for your help again !

 

[npoumaere]

Hi,

"calendar = 'thirty_day'" gives you a 360 days per year calendar with 12 months of 30 days each.
That's what I use as I am working on a highly idealized simulation, for which I'm not using exact daily forcing data taken from e.g. ERA5.
"NOLEAP" gives a 365 days per year and corresponds to our "real life" calendar with leap years removed.
"GREGORIAN" and "JULIAN" account for leap years (our "real life" calendar is gregorian), and are useful when you want to assimilate real life observations into your model.

Please have a look at e.g. MOM6 — DART 11.21.2 documentation or Time Ranges and Calendars | NA-CORDEX as I’m not very familiar with these specific questions myself.

Nelson

 

[ruyin Wang]

Hi,

"calendar = 'thirty_day'" gives you a 360 days per year calendar with 12 months of 30 days each.
That's what I use as I am working on a highly idealized simulation, for which I'm not using exact daily forcing data taken from e.g. ERA5.
"NOLEAP" gives a 365 days per year and corresponds to our "real life" calendar with leap years removed.
"GREGORIAN" and "JULIAN" account for leap years (our "real life" calendar is gregorian), and are useful when you want to assimilate real life observations into your model.

Please have a look at e.g. MOM6 — DART 11.21.2 documentation or Time Ranges and Calendars | NA-CORDEX as I’m not very familiar with these specific questions myself.

Nelson

Thank you very much for your help — I believe I now have a much better understanding.

I have a small follow-up question. I noticed that when using CORE climatological forcing, the calendar is set to “NOLEAP,” which allows the model to automatically cycle through a single year of forcing data. However, for JRA forcing, the default calendar is “Julian,” which seems to require the forcing time to strictly match the model time.

I was wondering whether your experiments involve this kind of “normal year forcing” setup. Specifically, in MOM/SIS, does the calendar type of the forcing files determine whether the model can loop over the forcing data?

Thank you again for your time and assistance.

Best regards,
Ruyin