[Mitgcm-support] Unstable 1deg Atl. model

[mitgcm-support at dev.mitgcm.org]

Hi. 

I'm building a 1deg resolution version of the 1/6deg North Atlantic model
(Helen supplied the original code). It's on the same rotated grid and uses
vector invariant momentum. KPP is on. I'm having difficulty finding a set
of dissipation parameters to give a stable code.

With the following data file I can't integrate more than a few thousand
steps (before the solver diverges). Increasing vertical viscosity or
bottom drag does not make much difference (decreasing them makes it much
less stable). Increasing horizontal dissipation makes it less stable.

Based on the verification examples my parameter settings look reasonable.
Any ideas why it's such a problem?  (Do I need CD?, is there something
special about vec inv mom?  Do I have to customize the bathymetry
somehow?).  I'm using NCEP reanalysis forcing which reaches a peak stress
of ~0.35Pa in winter, but that seems acceptable (incidentally, I rotated
the stress onto the rotated grid properly if anyone's interested). I'd
like to use GM eventually...

Thanks, Tom.

# ====================
# | Model parameters |
# ====================
#
# Continuous equation parameters
 &PARM01
 viscAh= 8.0E2,
 viscAz=1.0E-3,
 diffKhT= 8.0E2,
 diffKzT=3.0E-5,
 diffKhS= 8.0E2,
 diffKzS=3.0E-5,
 no_slip_sides=.FALSE.,
 eosType='POLY3',
 hFacMin=0.2,
 hFacMinDz=50.,
 cosPower=3.,
 staggerTimeStep=.TRUE.,
 vectorInvariantMomentum=.TRUE.,
 tempAdvScheme       = 4,
 saltAdvScheme       = 4,
 multiDimAdvection   = .FALSE.,
 implicitDiffusion   = .TRUE.,
 implicitViscosity   = .TRUE.
 bottomDragQuadratic = 0.001,
 &

# Elliptic solver parameters
 &PARM02
 cg2dMaxIters=1000,
 cg2dTargetResidual=1.E-13,
 &

# Time stepping parameters
 &PARM03
 nIter0  =0, 
 nTimeSteps= 34560, 
 deltaTmom=    1200.0, 
 deltaTClock = 1200.0,
 abEps=0.1,
 dumpFreq= 2592000.0, 
 movieFreq = 86400.0,
 tauThetaClimRelax=5184000.0, 
 tauSaltClimRelax=5184000.0, 
 &

# Gridding parameters
 &PARM04
 usingCartesianGrid=.FALSE.,
 usingSphericalPolarGrid=.TRUE.,
 delZ= 25.D0, 25.D0, 25.D0, 25.D0, 
       35.D0, 50.D0, 75.D0, 100.D0, 150.D0, 200.D0, 275.D0,
       350.D0, 415.D0, 450.D0, 500.D0, 500.D0, 500.D0, 500.D0,
       500.D0, 500.D0, 500.d0, 
 phiMin=-50,
 thetaMin=-74.75, 
 delY=112*1,
 delX=96*1,
 &

# Input datasets
 &PARM05
 hydrogThetaFile='theta.ic.bin.1deg',
 hydrogSaltFile ='salt.ic.bin.1deg',
 bathyFile      ='custom_bathy.bin.1deg',
 zonalWindFile  ='ncep.taus.1950-01-15.30day.bin',
 meridWindFile  ='ncep.taur.1950-01-15.30day.bin',
 surfQfile      ='ncep.net_hflx.1950-01-15.30day.bin',
 EmPmRfile      ='ncep.EmPmR.1950-01-15.30day.bin',
 thetaClimFile  ='SST.surf_relax.bin.1deg',
 saltClimFile   ='SSS.surf_relax.bin.1deg',
 fCoriogdegFile ='g_lat.bin.1deg',
 fCoriocdegFile ='c_lat.bin.1deg',
 &

# Open-boundaries
 &PARM06
 &