[MITgcm-devel] about dst3 (advection scheme 30)

[Martin Losch]

Hi there,

while searching for a possible bug that causes the noise in my 1/6th  
degree simulations, I came across a (probably historical) oddity in  
the DST3 routines gad_dst3_adv_x/y/r.F:
The original implementation was (only for the horizontal x-direction):
>         Rjp=(tracer(i+1,j)-tracer(i,j))*maskW(i+1,j,k,bi,bj)
>         Rj =(tracer(i,j)-tracer(i-1,j))*maskW(i,j,k,bi,bj)
>         Rjm=(tracer(i-1,j)-tracer(i-2,j))*maskW(i-1,j,k,bi,bj)
>
>         cfl=uVel(i,j,k,bi,bj)*deltaT*recip_dxc(i,j,bi,bj)
>         d0=(2.-abs(cfl))*(1.-abs(cfl))*oneSixth
>         d1=(1.-cfl)*(1.+cfl)*oneSixth
>         uT(i,j)=
>      &   0.5*(uTrans(i,j)+abs(uTrans(i,j)))
>      &      *( Tracer(i-1,j) + d0*Rj + d1*Rjm )
>      &  +0.5*(uTrans(i,j)-abs(uTrans(i,j)))
>      &      *( Tracer( i ,j) - d0*Rj + d1*Rjp )

Then Alistair changed it to make it look more like the flux limited  
version DST3FL (advection Scheme 33):
>         Rjp=(tracer(i+1,j)-tracer(i,j))*maskW(i+1,j,k,bi,bj)
>         Rj =(tracer(i,j)-tracer(i-1,j))*maskW(i,j,k,bi,bj)
>         Rjm=(tracer(i-1,j)-tracer(i-2,j))*maskW(i-1,j,k,bi,bj)
>
>         cfl=abs(uVel(i,j,k,bi,bj)*deltaT*recip_dxc(i,j,bi,bj))
>         d0=(2.-cfl)*(1.-cfl)*oneSixth
>         d1=(1.-cfl*cfl)*oneSixth
> c       thetaP=0.
> c       IF (Rj.NE.0.) thetaP=Rjm/Rj
>         thetaP=Rjm/(1.D-20+Rj)
>         psiP=d0+d1*thetaP
> c       psiP=max(0.,min(min(1.,psiP),(1.-cfl)/(1.D-20+cfl)*thetaP))
>         thetaM=Rjp/(1.D-20+Rj)
> c       thetaM=0.
> c       IF (Rj.NE.0.) thetaM=Rjp/Rj
>         psiM=d0+d1*thetaM
> c       psiM=max(0.,min(min(1.,psiM),(1.-cfl)/(1.D-20+cfl)*thetaM))
>         uT(i,j)=
>      &   0.5*(uTrans(i,j)+abs(uTrans(i,j)))
>      &      *( Tracer(i-1,j) + psiP*Rj )
>      &  +0.5*(uTrans(i,j)-abs(uTrans(i,j)))
>      &      *( Tracer( i ,j) - psiM*Rj )

The CVS comment is "Slight re-write to help debug flux limited form."
This version involves divisions by Rj, the case of Rj=0 is caught by  
adding a small number (1D-20), Patrick then changed this code again  
to make it adjoinable, but making the small number adjoinable. As the  
debugging of the flux limited form is probably finished, wouldn't it  
make more sense to revert the code to the orginal version which did  
not involve any division and should be adjoinable much more easily.  
Also, the original version is actually there, ifdef ALLOW_MATRIX. It  
all doesn't make too much sense to me. I would suggest reverting to  
the orgininal version, remove the ALLOW_MATRIX-ifdefs, and maybe add  
a flag such as ALLOW_DEBUGGING_DST3FL which contains the current  
version of the code. This modification should not affect any  
verification experiment, as the only one using advScheme=30 is  
matrix_example, where ALLOW_MATRIX is define, and even there it is  
"only" in the passive tracer package:
> csysm3::verification> grep 30, */input/data*
> matrix_example/input/data.ptracers: PTRACERS_advScheme(1)=30,

What do you think about this?

Martin