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

[Samar Khatiwala]

Martin, some history here...

I don't know about any "original version", but the code within ifdef  
ALLOW_MATRIX / endif was
something I wrote. I wanted a linear advection scheme (as DST3 should  
be), and the way the code
was previously implemented made it nonlinear. (A bit odd for a scheme  
that was NOT flux limiting.)
I got tired of using my own mods every time I wanted to use DST3, so  
Patrick kindly incorporated my
modifications in gad_dst3_adv_* within the ALLOW_MATRIX block. As  
currently implemented, the
(correct) linear version of DST3 is only executed when using the  
matrix package.

Samar

On Apr 2, 2006, at 7:56 AM, Martin Losch wrote:

> 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
> _______________________________________________
> MITgcm-devel mailing list
> MITgcm-devel at mitgcm.org
> http://mitgcm.org/mailman/listinfo/mitgcm-devel