--- MITgcm/pkg/generic_advdiff/gad_dst3_adv_y.F	2004/09/24 16:53:46	1.4
+++ MITgcm/pkg/generic_advdiff/gad_dst3_adv_y.F	2005/04/03 16:05:34	1.5
@@ -1,10 +1,10 @@
-C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/generic_advdiff/gad_dst3_adv_y.F,v 1.4 2004/09/24 16:53:46 jmc Exp $
+C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/generic_advdiff/gad_dst3_adv_y.F,v 1.5 2005/04/03 16:05:34 heimbach Exp $
 C $Name:  $
 
 #include "GAD_OPTIONS.h"
 
       SUBROUTINE GAD_DST3_ADV_Y( 
-     I           bi,bj,k,deltaT,
+     I           bi,bj,k,deltaTloc,
      I           vTrans, vVel,
      I           maskLocS, tracer,
      O           vT,
@@ -19,11 +19,13 @@
 C     == GLobal variables ==
 #include "SIZE.h"
 #include "GRID.h"
+#include "EEPARAMS.h"
+#include "PARAMS.h"
 #include "GAD.h"
 
 C     == Routine arguments ==
       INTEGER bi,bj,k
-      _RL deltaT
+      _RL deltaTloc
       _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
       _RL vVel(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
       _RS maskLocS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
@@ -37,6 +39,14 @@
       _RL Rjm,Rj,Rjp,cfl,d0,d1
       _RL psiP,psiM,thetaP,thetaM
       _RL vFld
+      _RL smallNo
+      _RL Rjjm,Rjjp
+
+      IF (inAdMode) THEN
+       smallNo = 1.0D-20
+      ELSE
+       smallNo = 1.0D-20
+      ENDIF
 
       DO i=1-Olx,sNx+Olx
        vT(i,1-Oly)=0.
@@ -52,19 +62,25 @@
 c       vFld = vVel(i,j,k,bi,bj)
         vFld = vTrans(i,j)*recip_dxG(i,j,bi,bj)
      &       *recip_drF(k)*recip_hFacS(i,j,k,bi,bj)
-        cfl=abs(vFld*deltaT*recip_dyC(i,j,bi,bj))
+        cfl=abs(vFld*deltaTloc*recip_dyC(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))
+        IF ( ABS(Rj).LT.smallNo .OR.
+     &       ABS(Rjm).LT.smallNo ) THEN
+         thetaP=0.
+         psiP=0.
+        ELSE   
+         thetaP=(Rjm+smallNo)/(smallNo+Rj)
+         psiP=d0+d1*thetaP
+        ENDIF
+        IF ( ABS(Rj).LT.smallNo .OR.
+     &       ABS(Rjp).LT.smallNo ) THEN
+         thetaM=0.
+         psiM=0.
+        ELSE
+         thetaM=(Rjp+smallNo)/(smallNo+Rj)
+         psiM=d0+d1*thetaM
+        ENDIF
         vT(i,j)=
      &   0.5*(vTrans(i,j)+abs(vTrans(i,j)))
      &      *( Tracer(i,j-1) + psiP*Rj )