[MITgcm-devel] viscosity on CS grid

[Jean-Michel Campin]

Hi,

it seems to me that we have 2 problems with horizontal 
viscosity on the cubed-sphere grid, at the corners:

1) Leith scheme (C2 or C4):
for example, SW corner:
vort3(0,1) is not computed (set to zero).
but is used to compute viscAh_D(0,1) and viscAh_Z(1,1)
and later uDissisp(1,1).
This means that we will get different values
on the 2 faces for the same physical point on the edge 
of the cube, near a corner.

2) biharmonic viscosity:
generally requires Olx=3,Oly=3.
but around the corner, we have some problems to compute
del2u,del2v on a wide enough stencil. At the end
uDissisp & vDissisp are not correct (different values
on the 2 faces for the same physical point on the edge 
of the cube, near a corner).

I can see 2 type of solutions:
a) add an exchange of del2u & del2v to fix the biharmonic
problem. Unfortunately, del2u &  del2v are only 2D local
arrays, inside bi,bj loop. It means that the need to split
the bi,bj loops, may be also the k loop, a pain.
To fix the Leith scheme, we need to exchange the vorticity (vort3),
same problem with local array inside bi,bj loops, but in addition,
we still don't have the routine to do the exchange at a corner
of the C-grid.
Apart from the modifications needed in the code, those additional
exchanges will slow down the model.

b) try to fix the Leith scheme without additional exchange,
by computing the relative vorticity over a wider stencil.
I started to change EXCH_UV and it doesn't seem impossible
to make it work (still few problems with grid length
and grid cell area).
Now, if this works, this will reduce the biharmonic problem
to only the gradient of divergence part of the flow (since the 
vorticity will be available where it's needed), and this can be 
fixed without additional exchanges, like the multidimension advection
(for each direction, do a local copy into the corner-halo region 
from the same tile, before computing each gradient component).

Comments/suggestions ?

Jean-Michel