[MITgcm-devel] ggl90_calc

[Martin Losch]

Hi Jean-Michel,

I am confused by my own code. I need to revise a few of my recent  
comments:
On Oct 31, 2009, at 4:31 AM, Jean-Michel Campin wrote:

> Hi Martin and David,
>
> I look at S/R ggl90_calc.F (more precisely, to the implicit diffusion
> of TKE) and noticed few thing that I don't like too much
> (at least misleading, but some might also be more like a bugg).
>
> It's all checked-in as comments in the code ("C- jmc:"):
>> danton{ggl90}% egrep '^C- jmc:' ggl90_calc.F
>> C- jmc: concerned about missing a deltaT since gTKE is already the  
>> future TKE.
This is still a bug
>> C- jmc: concerned that a(k=2) should always be zero
I do not get this comment. a(k=1) = 0. This is the surface value and  
is moved to the rhs (gTKE(k=1)), why do you want to move the surface  
value to k=2? why should a(k=2) = 0 and likewise why should gTKE(k=2)  
be updated?
>> C- jmc: concerned that c(k) from k=kLow to k=Nr should always be zero
I thought that this is achieved through maskC(k)*maskC(k-1)
>> C- jmc: this is dangerous since klowC could be zero if land column
True, what we really need is this
           kp1=k+1
           if ( klowC(i,j,bi,bj) .gt. 0 ) kp1=min(klowC(i,j,bi,bj),kp1)
>> C- jmc: concerned about conservation when a or c are changed after  
>> computing b
this can easily be changed by moving the boundary conditions loop up  
before the compuation of b
>> C- jmc: would be much better to update the provisional TKE (i.e.  
>> gTKE) at k=2
see above.

So in summary I see only too problems (please point me to the things  
that I do not see, most likely the misleading part):
1. bug in the time stepping of the horizontal diffusivities
2. bug in the computation of kp1 for klowc = 0

Then ggl90_solve.F is really a copy of solve_tridiagonal (with small  
difference) and could be removed/replaced. Unfortunately switching  
between the routines cause an error (8 digits) and I am not clear if  
this is only numerical truncation.

Martin