Checkerboard instability

[dlebars]

Dear all,

I should be doing something wrong but I can't find out so any help is welcome.
I develop a 3D ROMS (version 3.7) configuration of the european coast forced with ERA-interim fields and bulk formula with horizontal resolution of 20km and 30 vertical layers. I get what looks like a checkerboard instability growing within a few days leading the model to blow up. I tried different boundary conditions, including closing all boundaries and this has no impact on the instability. I tried different viscosities, harmonic, biharmonic, smagorinsky the only way to get rid of the instability is to put the viscosity around 100 times bigger than what I would except at this resolution (for instance a harmonic viscosity of 10^5 m^2/s is necessary...). I also tried smoothing the bottom topography a lot with no impact.

Please find bellow an example plot of the instability and part of the log file.

Thank you for your help,
Dewi

### parts of the log file ###

Physical Parameters, Grid: 01
=============================

291840 ntimes Number of timesteps for 3-D equations.
30.000 dt Timestep size (s) for 3-D equations.
3 ndtfast Number of timesteps for 2-D equations between
each 3D timestep.
1 ERstr Starting ensemble/perturbation run number.
1 ERend Ending ensemble/perturbation run number.
0 nrrec Number of restart records to read from disk.
T LcycleRST Switch to recycle time-records in restart file.
28800 nRST Number of timesteps between the writing of data
into restart fields.
1 ninfo Number of timesteps between print of information
to standard output.
T ldefout Switch to create a new output NetCDF file(s).
720 nHIS Number of timesteps between the writing fields
into history file.
1 ntsAVG Starting timestep for the accumulation of output
time-averaged data.
28800 nAVG Number of timesteps between the writing of
time-averaged data into averages file.
1 ntsDIA Starting timestep for the accumulation of output
time-averaged diagnostics data.
720 nDIA Number of timesteps between the writing of
time-averaged data into diagnostics file.
1.0000E+02 nl_tnu2(01) NLM Horizontal, harmonic mixing coefficient
(m2/s) for tracer 01: temp
1.0000E+02 nl_tnu2(02) NLM Horizontal, harmonic mixing coefficient
(m2/s) for tracer 02: salt
1.0000E+03 nl_visc2 NLM Horizontal, harmonic mixing coefficient
(m2/s) for momentum.
F LuvSponge Turning OFF sponge on horizontal momentum.
F LtracerSponge(01) Turning OFF sponge on tracer 01: temp
F LtracerSponge(02) Turning OFF sponge on tracer 02: salt
1.0000E-05 Akt_bak(01) Background vertical mixing coefficient (m2/s)
for tracer 01: temp
1.0000E-05 Akt_bak(02) Background vertical mixing coefficient (m2/s)
for tracer 02: salt
1.0000E-04 Akv_bak Background vertical mixing coefficient (m2/s)
for momentum.
3.0000E-04 rdrg Linear bottom drag coefficient (m/s).
3.0000E-03 rdrg2 Quadratic bottom drag coefficient.
2.0000E-02 Zob Bottom roughness (m).
2.0000E+00 blk_ZQ Height (m) of surface air humidity measurement.
2.0000E+00 blk_ZT Height (m) of surface air temperature measurement.
1.0000E+01 blk_ZW Height (m) of surface winds measurement.
5 lmd_Jwt Jerlov water type.
2 Vtransform S-coordinate transformation equation.
4 Vstretching S-coordinate stretching function.
7.0000E+00 theta_s S-coordinate surface control parameter.
2.0000E+00 theta_b S-coordinate bottom control parameter.
50.000 Tcline S-coordinate surface/bottom layer width (m) used
in vertical coordinate stretching.
1025.000 rho0 Mean density (kg/m3) for Boussinesq approximation.
16451.000 dstart Time-stamp assigned to model initialization (days).
19480101.00 time_ref Reference time for units attribute (yyyymmdd.dd)
3.6000E+02 Tnudg(01) Nudging/relaxation time scale (days)
for tracer 01: temp
3.6000E+02 Tnudg(02) Nudging/relaxation time scale (days)
for tracer 02: salt
3.6000E+02 Znudg Nudging/relaxation time scale (days)
for free-surface.
3.6000E+02 M2nudg Nudging/relaxation time scale (days)
for 2D momentum.
3.6000E+02 M3nudg Nudging/relaxation time scale (days)
for 3D momentum.
3.6000E+01 obcfac Factor between passive and active
open boundary conditions.
T VolCons(1) NLM western edge boundary volume conservation.
T VolCons(2) NLM southern edge boundary volume conservation.
T VolCons(3) NLM eastern edge boundary volume conservation.
T VolCons(4) NLM northern edge boundary volume conservation.
10.000 T0 Background potential temperature (C) constant.
35.000 S0 Background salinity (PSU) constant.
-1.000 gamma2 Slipperiness variable: free-slip (1.0) or
no-slip (-1.0).

Lateral Boundary Conditions: NLM
============================

Variable Grid West Edge South Edge East Edge North Edge
--------- ---- ---------- ---------- ---------- ----------

zeta 1 Closed Closed Closed Closed

ubar 1 Closed Closed Closed Closed

vbar 1 Closed Closed Closed Closed

u 1 Closed Closed Closed Closed

v 1 Closed Closed Closed Closed

temp 1 Closed Closed Closed Closed

salt 1 Closed Closed Closed Closed

Activated C-preprocessing Options:

NORTH_SEA4 Extended North Sea, 0.25deg Resolution
ANA_BSFLUX Analytical kinematic bottom salinity flux.
ANA_BTFLUX Analytical kinematic bottom temperature flux.
ANA_SRFLUX Analytical kinematic shortwave radiation flux.
ASSUMED_SHAPE Using assumed-shape arrays.
AVERAGES Writing out time-averaged nonlinear model fields.
BULK_FLUXES Surface bulk fluxes parameterization.
CURVGRID Orthogonal curvilinear grid.
DIAGNOSTICS_UV Computing and writing momentum diagnostic terms.
DIURNAL_SRFLUX Modulate shortwave radiation by the local diurnal cycle.
DOUBLE_PRECISION Double precision arithmetic.
EMINUSP Compute Salt Flux using E-P.
LMD_BKPP KPP bottom boundary layer mixing.
LMD_CONVEC LMD convective mixing due to shear instability.
LMD_DDMIX LMD double-diffusive mixing.
LMD_MIXING Large/McWilliams/Doney interior mixing.
LMD_NONLOCAL LMD convective nonlocal transport.
LMD_RIMIX LMD diffusivity due to shear instability.
LMD_SKPP KPP surface boundary layer mixing.
LONGWAVE Compute net longwave radiation internally.
MASKING Land/Sea masking.
MIX_ISO_TS Mixing of tracers along isopycnal surfaces.
MIX_S_UV Mixing of momentum along constant S-surfaces.
MPI MPI distributed-memory configuration.
NONLINEAR Nonlinear Model.
NONLIN_EOS Nonlinear Equation of State for seawater.
POWER_LAW Power-law shape time-averaging barotropic filter.
PRSGRD31 Standard density Jacobian formulation (Song, 1998).
PROFILE Time profiling activated .
RI_HORAVG Smooth Richardson number horizontally.
RI_VERAVG Smooth Richardson number vertically.
RHO_SURF Include difference between rho0 and surface density.
!RST_SINGLE Double precision fields in restart NetCDF file.
SALINITY Using salinity.
SOLVE3D Solving 3D Primitive Equations.
TS_C4HADVECTION Fourth-order centered horizontal advection of tracers.
TS_C4VADVECTION Fourth-order centered vertical advection of tracers.
TS_DIF2 Harmonic mixing of tracers.
UV_ADV Advection of momentum.
UV_COR Coriolis term.
UV_U3HADVECTION Third-order upstream horizontal advection of 3D momentum.
UV_C4VADVECTION Fourth-order centered vertical advection of momentum.
UV_QDRAG Quadratic bottom stress.
UV_VIS2 Harmonic mixing of momentum.
VAR_RHO_2D Variable density barotropic mode.

[jcwarner]

i cant download the attachment because my computer is being overly secure. but did you activate
#define N2S2_HORAVG

-j

[dlebars]

Dear jcwarner,

For the moment I am using LMD_MIXING so from what I understand N2S2_HORAVG is not used in this case.
But I am going to try GLS_MIXING with N2S2_HORAVG maybe it solves my issue.

Thanks for the tip,
Dewi

[dlebars]

Using GLS vertical mixing parameterization with standard k-epsilon implementation and N2S2_HORAVG does not remove the instability.

[kate]

It seems to me that for 20 km, you ought to be able to use a longer timestep. This is what I use for an 11 km domain:

Code: Select all

          DT == 225.0d0
     NDTFAST == 20

also

Code: Select all

        TNU2 == 5.0d0  5.0d0                    ! m2/s
       VISC2 == 25.0d0                          ! m2/s

plus using the third-order upwind advection to keep things smoother:

Code: Select all

# define TS_U3HADVECTION
# define TS_C4VADVECTION

ETA: Woah, bizarre plot. Can you plot it earlier to see where it originates from? It appears to be in all the deep water, so I'd try a shorter barotropic timestep.

[jklinck]

A checkerboard instability pattern is due horizontal diffusion/viscosity. Look at your time step relative to the grid spacing and horizontal diffusivity. With laplacian horizontal diffusivity A dt / dx^2 needs to be less than 1 (maybe 1/2) depending on the time stepping scheme. With variable horizontal diffusivity choices, look for the largest diffusivity/viscosity. Vertical diffusion is done implicitly, so it is unlikely that this pattern is due to a vertical subgridscale model.

[wilkin]

3 ndtfast Number of timesteps for 2-D equations between

You only have 3 time steps of the barotropic mode for every baroclinic time step. So the S-shaped time filter that averages fast-time-step (barotropic) zeta to the slow time step has very few weights. Most users have NDTFAST no less than 20, and typically 30 or more.

PRSGRD31 Standard density Jacobian formulation (Song, 1998).

The vast majority of the test cases use the splines density Jacobian (Shchepetkin, 2000) activated with #define DJ_GRADPS. I suggest you switch to this option to decrease pressure gradient truncation errors.

But I'm pretty sure your problem is NDTFAST.

[dlebars]

The problem is solved! It was indeed NDTFAST, I sought it was the barotropic time step in seconds so I was decreasing it to remove the instability which made things worse.
Thank you for your help!