[Mitgcm-support] Re: upwind

[mitgcm-support at dev.mitgcm.org]

Alistair,

the 3. order upwind scheme is still not running but crashes at the
northern boundary. Even the first order has its problems and is very
noisy - this seems to be due to the lopped cells options (the major
difference to the global model where first order upwind was fine). But
compared to 3. order upwind the first order scheme does not crash and is
noisy only in the level where the daily wind is applied.

Could you please have a look at some results at
http://www.ecco.ucsd.edu/internal/subpol/comparison.html
and write some suggestions?

I might be in a complete new parameter range than before, and the lopped
cells option might even narrow it. I am currently running a test without
lopped cells to isolate the problem.

Thanks a lot.
Ciao,
	Arne

Alistair Adcroft wrote:
> 
> Arne Biastoch wrote:
> > (i) how do you get the numbers -0.125 and 0.5? When I look at the 3.
> > order
> > upwind discretization (e.g. 2.28 in Haidvogel and Beckmann), there are
> > factors of -1/3 and 4/3. Even if you consider the divisions by 2dx and
> > 4dx
> > the remaining factors are -1/12 and 2/3.
> 
> Assuming U>0 then the stencil for the right hand flux is
> composed of centered average (4/8=1/2=0.5)
>    [  0    0   4/8  4/8 ]
> and a side ways diffusion (1/8=0.125)
>  + [  0  -1/8  2/8 -1/8 ]
> which gives a total stencil of
>  = [ 0   -1/8  6/8  3/8 ]  = R
> Similarly, the left hand flux has a stencil of
>    [  0   4/8  4/8   0  ]
>  + [-1/8  2/8 -1/8   0  ]
>  = [-1/8  6/8  3/8   0  ]  = L
> If we now take the difference between the fluxes (div.F)
> and add up the stencils we get
> R-L =
>    [1/8  -7/8  3/8  3/8 ]
> which is what is written out in 2.27 in Haidvogel and Beckmann.
> 
> Eq. 2.28 is fourth order which can be written as a
> centered flux - a centered evaluation of d_xx times 1/6.
>   R = [  0     0   6/12  6/12 0 ] + [   0   -1/12 1/12  1/12 -1/12 ]
>   L = [  0    6/12 6/12   0   0 ] + [ -1/12  1/12 1/12 -1/12   0   ]
> R-L = [ 1/12 -8/12  0    8/12 -1/12 ]
> We write things like this so we can evaluate boundary conditions
> more naturally.
> 
> > (ii) for the case uTrans < 0 TxM is always zero.
> >            TxM=(theta(i-1,j,k,bi,bj)-theta(i-1,j,k,bi,bj))*
> >      &          _maskW(i,j,k,bi,bj)
> 
> bug fix:
>              TxM=(theta(i,j,k,bi,bj)-theta(i-1,j,k,bi,bj))*
>        &          _maskW(i,j,k,bi,bj)
> 
> Alistair.
> cc support

-- 

Dr. Arne Biastoch
Scripps Institution of Oceanography        phone: +1-858-822-3787
Physical Oceanography Research Division    fax  : +1-858-534-9820
MS: 0230                                   email: abiastoch at ucsd.edu
8605 La Jolla Shores Dr.            
La Jolla, CA 92093-0230, U.S.A.   http://www.ecco.ucsd.edu/~biastoch