--- MITgcm/pkg/generic_advdiff/gad_dst3fl_adv_x.F	2001/09/05 17:46:03	1.2
+++ MITgcm/pkg/generic_advdiff/gad_dst3fl_adv_x.F	2001/09/10 00:24:14	1.3
@@ -1,4 +1,4 @@
-C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/generic_advdiff/gad_dst3fl_adv_x.F,v 1.2 2001/09/05 17:46:03 heimbach Exp $
+C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/generic_advdiff/gad_dst3fl_adv_x.F,v 1.3 2001/09/10 00:24:14 adcroft Exp $
 C $Name:  $
 
 #include "GAD_OPTIONS.h"
@@ -32,7 +32,7 @@
 
 C     == Local variables ==
       INTEGER i,j
-      _RL Rjm,Rj,Rjp,cfl,d0,d1,psiP,thetaP,psiM,thetaM
+      _RL Rjm,Rj,Rjp,cfl,d0,d1,psiP,psiM,thetaP,thetaM
 
       DO j=1-Oly,sNy+Oly
        uT(1-Olx,j)=0.
@@ -43,26 +43,24 @@
         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
-        thetaP=Rjm/(1.D-30+Rj)
+        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
-        psiP=max(0., min(min(1.,psiP),(1.-cfl)/cfl*thetaP) )
-        thetaM=Rjp/(-1.D-30+Rj)
+        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
-        psiM=max(0., min(min(1.,psiM),(1.-cfl)/cfl*thetaM) )
+        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 )
-cph *note* put these comments after end of continued line
-cph        to ensure TAMC compatibility
-c    &   0.5*(uTrans(i,j)+abs(uTrans(i,j)))
-c    &      *( Tracer(i-1,j) + d0*Rj + d1*Rjm )
-c    &  +0.5*(uTrans(i,j)-abs(uTrans(i,j)))
-c    &      *( Tracer( i ,j) - d0*Rj + d1*Rjp )
 
        ENDDO
       ENDDO