[MITgcm-support] Re: viscosity questions

Mon Jun 26 10:34:49 EDT 2006

Hi Gianmaria (and mitgcm_support for others who may be interested in  
how to set viscosities initially),
   The viscosity code has been very stable for the hydrostatic part  
of the code for a while now.  Unless you are intending to run the  
nonhydrostatic version (which will probably not do anything at 1/8th  
degree), then you should be OK with the month-old version.  I am  
pretty sure we had all the bugs out by checkpoint57x or later.  If  
you want to use the nonhydrostatic version, you will need to upgrade  
to the most recent CVS check-in.
   I have a few suggestions as to how to get started.  If you want to  
use KPP and GMREDI for diffusion and minimize across-isopycnal  
diffusion, I'd use the 30 advection scheme or the flux-limited 33  
advection scheme and eliminate all diffKhT and diffK4T.  That should  
be stable, and while there will be cross-ispycnal implicit mixing due  
to the advection scheme, it should be as small as possible.

With these schemes, you should be careful to turn on the diagnostics:
   86 |ADVx_TH | 40 |UU   087MR      |degC.m^3/s      |Zonal       
Advective Flux of Pot.Temperature
   87 |ADVy_TH | 40 |VV   086MR      |degC.m^3/s      |Meridional  
Advective Flux of Pot.Temperature

   (if you're using salt)
   93 |ADVx_SLT| 40 |UU   094MR      |psu.m^3/s       |Zonal       
Advective Flux of Salinity
   94 |ADVy_SLT| 40 |VV   093MR      |psu.m^3/s       |Meridional  
Advective Flux of Salinity

   89 |DIFx_TH | 40 |UU   090MR      |degC.m^3/s      |Zonal       
Diffusive Flux of Pot.Temperature
   90 |DIFy_TH | 40 |VV   089MR      |degC.m^3/s      |Meridional  
Diffusive Flux of Pot.Temperature

   As far as viscosity, you should probably be using a Leith  
viscosity at 1/8th degree, and you can experiment with whether using  
LeithD is an improvement (often it reduces noise in W and may make  
your time-steps able to be bigger).  You could also try Leith 
+Smagorinsky which is a bit more theoretically satisfying, but may  
not have the full numerical advantages of LeithD.  You may want to  
turn on the diagnostics:
  101 |VISCAHD | 40 |SM      MR      |m^2/s           |Harmonic  
Viscosity Coefficient (m2/s) (Div Pt)
Which will let you know what value of viscosity was being used.

Here'e the basic way of setting up the model.  These are testing  
runs, so you don't have to get everything else right to get these  
right.  Also, you don't have to spin up too long, since the  
barotropic flow will spin up quickly O(a year), and getting the rough  
tuning right should be ok without a spun-up baroclinic flow.

Setup #1:
Leith Only:
1) Set vishAhGridMax=1 (or maybe 0.5 to be safe, no real difference  
there.)
2) Set viscC2Leith=1.
3) All other viscosities will default to zero (or you could set them  
to a tiny value if you like).
4) Do a trial run and look at the output.  Pay particular attention  
to the boundary currents.  Is their width acceptable?  You may want  
to look at the diagnostics:
   132 |Um_Diss | 40 |UU   133MR      |m/s^2           |U momentum  
tendency from Dissipation
   133 |Vm_Diss | 40 |VV   132MR      |m/s^2           |V momentum  
tendency from Dissipation
   134 |Um_Advec| 40 |UU   135MR      |m/s^2           |U momentum  
tendency from Advection terms
   135 |Vm_Advec| 40 |VV   134MR      |m/s^2           |V momentum  
tendency from Advection terms
In theory, the advective part should dominate over most of the  
domain, but in a few regions the dissipation terms may be important.   
Also, look for regions where the field is checkerboard patterned.   
You may want to look at W (which will be the noisiest).
5) Change viscC2Leith up or down a bit to improve. Go to 4)

Setup #2:
Leith+LeithD (May be necessary if checkerboard noise in W doesn't go  
away with Leith above, may also make larger timesteps possible).
1) Set vishAhGridMax=1 (or maybe 0.5 to be safe, no real difference  
there.)
2) Set viscC2Leith=viscC2LeithD=1.
3) All other viscosities will default to zero (or you could set them  
to a tiny value if you like).
4) Do a trial run and look at the output.  You may want to check out  
the viscosities:
   113 |VAHDLTH | 40 |SM      MR      |m^2/s           |Leith Harm  
Visc Coefficient (m2/s) (Div Pt)
   117 |VAHDLTHD| 40 |SM      MR      |m^2/s           |LeithD Harm  
Visc Coefficient (m2/s) (Div Pt)
5) Change viscC2Leith and viscC2LeithD keeping them equal. Go to 4)

Setup #3:
Leith+Smagorinsky (Perhaps the most 'physical', but maybe overly  
complicated or computer-time consuming)
1) Set vishAhGridMax=1 (or maybe 0.5 to be safe, no real difference  
there.)
2) Set viscC2Leith=viscC2Smag=1.
3) All other viscosities will default to zero (or you could set them  
to a tiny value if you like).
4) Do a trial run and look at the output.  You may want to focus on  
the ratio of the viscosities from the diagnostic:
    113 |VAHDLTH | 40 |SM      MR      |m^2/s           |Leith Harm  
Visc Coefficient (m2/s) (Div Pt)
    121 |VAHDSMAG| 40 |SM      MR      |m^2/s           |Smagorinsky  
Harm Visc Coefficient (m2/s) (Div Pt)
Is either one doing most of the work?  Are regions that appear overly  
viscous dominated by one or the other?
5) Change viscC2Leith and viscC2Smag, but there is no reason to keep  
them equal. Go to 4)

Other considerations.
Biharmonic Viscosities:  You may prefer biharmonic viscosities (I  
personally don't).  Just do the same as above, except with  
ViscA4Gridmax and viscC4leith, etc.  If you are unable to eliminate  
checkerboard patterns with algorithms above, you can use biharmonic  
and harmonic simultaneously, e.g. viscC2leith=1, viscC4leith=0.1.   
This should add a bit of biharmonic at the gridscale.

Background constant:  You can also set a background constant  
viscosity (viscAh) or one that scales with resolution (viscAhgrid).   
These will get added to the above.  This might  be useful if you want  
to ensure that a particular boundary layer width (e.g., L_Munk= 
(viscAh/beta)^(1/3)) is resolved (approximately L_Munk=2*deltaX).

Boundary Conditions:  You can use either no_slip_sides=.true.  
or .false. with all of the above viscosities,  This will depend on  
you simulation, and what the results look like.  There are very  
important physical considerations as well (see, e.g., Fox-Kemper JPO  
2005).

Vertical Viscosities;  Currently, the vertical viscosities are  
independent of the considerations above.  Someday I hope to have a  
Smagorinsky-scaled vertical viscosity, but it's not done yet.  So,  
set you viscAr=1e-4 m^2/s or your favorite value, and carry on.  You  
will also need to decide on whether you will have no_slip_bottom and  
additional bottomDragLinear  or bottomDragQuadratic.  The latter can  
be useful in limiting the inverse cascade in models with flat bottoms.

I think that's pretty much it.
    Cheers,
      -Baylor