Temp/salinity and thickness bug when initializing with an ice shelf

[claireyung]

Hello, I've come across strange features in the initialization of temperature/salinity and layer thickness/interface levels when an ice shelf is present. I initialized the model with a simple triangular ice shelf and TS_CONFIG = "ISOMIP" and THICKNESS_CONFIG = "ISOMIP" (though the same happens with other config options, e.g. TS=linear or Thickness=uniform) which is meant to produce a linear stratification according to given temperature and salinity bounds. However, although it looks mostly linear, there are strange periodic deviations from linearity in the initial conditions file which are highlighted below in this figure, showing the anomaly of the salt initial conditions from the theoretical linear profile.

Sea level/layer thicknesses also have very slight perturbations from their expected smooth distributions (for sigma_shelf_zstar coordinates) when TRIM_IC_FOR_P_SURF is not used (turning it on seems to fix those discontinuities but not the T/S ones). When I run the model, horizontal velocities emerge very quickly with the same spatial frequency as the salt anomalies suggesting that the anomalous horizontal density gradients drive this unwanted flow.

A temporary fix: I can get around this issue by using netcdf files as input for the thicknesses (with ice shelf accounted for), and by asking the model to skip the calc_sfc_displacement step in MOM_state_initialization.F90. This somehow results in a linear T/S profile and very low (order 10^(-8) m/s) velocities for the simple triangular ice shelf, however requires running the model twice. If it is helpful, I have linked a github repo that includes the files/MOM_input/override used in my test cases.

[claireyung]

Verifying that this issue was resolved by 88d9eb7 and 25989ad

With the current dev/gfdl branch, the following parameters work well for me for ALE ice shelf initialisation:

#override REMAP_BOUNDARY_EXTRAP = False ! since we don't want remap extrapolation when the model is running
#override INIT_BOUNDARY_EXTRAP = True ! not needed for ALE z coordinates, but doesn't hurt to have it on
#override TRIMMING_USES_REMAPPING = True 
#override TRIM_IC_Z_TOLERANCE = 1.e-14 !

(may need to increase this tolerance depending on the geometry if model crashes (#609), but want it to be as small as possible)

Additionally, I use:

MIN_THICKNESS = 1E-12 !Small minimum thickness so viscosity acts to arrest flows in these layers
ICE_SHELF = True !Turn on ice shelf
ICE_PROFILE_CONFIG = “FILE” !Use input file of ice draft
ICE_THICKNESS_FILE = “XX.nc” ! with ICE_THICKNESS_VARNAME as the draft and ICE_AREA_VARNAME as the area
TRIM_IC_FOR_P_SURF = True !Use trim for ice submodule to add ice shelf
SURFACE_PRESSURE_FILE = “XX.nc” !You can use the same file as for the ice shelf draft
SURFACE_PRESSURE_VAR !Whatever ICE THICKNESS VARNAME is
SURFACE_PRESSURE_SCALE = 8820.0 !Or whatever G_EARTH × DENSITY_ICE is, so that you get pressure from the ice draft file XX.nc

Folder with some plots and configs showing the initialisation verification here if useful to someone in the future!