--- MITgcm/pkg/mom_fluxform/mom_fluxform.F	2005/09/23 15:19:38	1.26
+++ MITgcm/pkg/mom_fluxform/mom_fluxform.F	2005/11/24 00:06:38	1.30
@@ -1,4 +1,4 @@
-C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/mom_fluxform/mom_fluxform.F,v 1.26 2005/09/23 15:19:38 jmc Exp $
+C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/mom_fluxform/mom_fluxform.F,v 1.30 2005/11/24 00:06:38 heimbach Exp $
 C $Name:  $
 
 CBOI
@@ -145,6 +145,7 @@
       _RL  ArDvdrFac
       _RL  fvFac
       _RL  mtFacV
+      _RL  sideMaskFac
       LOGICAL bottomDragTerms,harmonic,biharmonic,useVariableViscosity
 CEOP
 
@@ -165,6 +166,11 @@
         tension(i,j)= 0.
         guDiss(i,j) = 0.
         gvDiss(i,j) = 0.
+#ifdef ALLOW_AUTODIFF_TAMC
+        vort3(i,j)   = 0. _d 0
+        strain(i,j)  = 0. _d 0
+        tension(i,j) = 0. _d 0
+#endif
        ENDDO
       ENDDO
 
@@ -193,6 +199,14 @@
         ArDvdrFac  = 0.
       ENDIF
 
+C note: using standard stencil (no mask) results in under-estimating
+C       vorticity at a no-slip boundary by a factor of 2 = sideDragFactor
+      IF ( no_slip_sides ) THEN
+        sideMaskFac = sideDragFactor
+      ELSE
+        sideMaskFac = 0. _d 0
+      ENDIF
+
       IF (     no_slip_bottom
      &    .OR. bottomDragQuadratic.NE.0.
      &    .OR. bottomDragLinear.NE.0.) THEN
@@ -232,10 +246,28 @@
       ENDDO
 
       CALL MOM_CALC_KE(bi,bj,k,2,uFld,vFld,KE,myThid)
-      CALL MOM_CALC_HDIV(bi,bj,k,2,uFld,vFld,hDiv,myThid)
-      CALL MOM_CALC_RELVORT3(bi,bj,k,uFld,vFld,hFacZ,vort3,myThid)
-      CALL MOM_CALC_TENSION(bi,bj,k,uFld,vFld,tension,myThid)
-      CALL MOM_CALC_STRAIN(bi,bj,k,uFld,vFld,hFacZ,strain,myThid)
+      IF ( momViscosity) THEN
+        CALL MOM_CALC_HDIV(bi,bj,k,2,uFld,vFld,hDiv,myThid)
+        CALL MOM_CALC_RELVORT3(bi,bj,k,uFld,vFld,hFacZ,vort3,myThid)
+        CALL MOM_CALC_TENSION(bi,bj,k,uFld,vFld,tension,myThid)
+        CALL MOM_CALC_STRAIN(bi,bj,k,uFld,vFld,hFacZ,strain,myThid)
+        DO j=1-Oly,sNy+Oly
+         DO i=1-Olx,sNx+Olx
+           IF ( hFacZ(i,j).EQ.0. ) THEN
+             vort3(i,j)  = sideMaskFac*vort3(i,j)
+             strain(i,j) = sideMaskFac*strain(i,j)
+           ENDIF
+         ENDDO
+        ENDDO
+#ifdef ALLOW_DIAGNOSTICS
+        IF ( useDiagnostics ) THEN
+          CALL DIAGNOSTICS_FILL(hDiv,   'momHDiv ',k,1,2,bi,bj,myThid)
+          CALL DIAGNOSTICS_FILL(vort3,  'momVort3',k,1,2,bi,bj,myThid)
+          CALL DIAGNOSTICS_FILL(tension,'Tension ',k,1,2,bi,bj,myThid)
+          CALL DIAGNOSTICS_FILL(strain, 'Strain  ',k,1,2,bi,bj,myThid)
+        ENDIF
+#endif
+      ENDIF
 
 C---  First call (k=1): compute vertical adv. flux fVerU(kUp) & fVerV(kUp)
       IF (momAdvection.AND.k.EQ.1) THEN
@@ -361,11 +393,11 @@
 
 C     Laplacian and bi-harmonic terms, Zonal  Fluxes -> fZon
         CALL MOM_U_XVISCFLUX(bi,bj,k,uFld,v4F,fZon,
-     I    viscAh_Z,viscAh_D,viscA4_Z,viscA4_D,myThid)
+     I    viscAh_D,viscA4_D,myThid)
 
 C     Laplacian and bi-harmonic termis, Merid Fluxes -> fMer
         CALL MOM_U_YVISCFLUX(bi,bj,k,uFld,v4F,hFacZ,fMer,
-     I    viscAh_Z,viscAh_D,viscA4_Z,viscA4_D,myThid)
+     I    viscAh_Z,viscA4_Z,myThid)
 
 C     Eddy component of vertical flux (interior component only) -> fVrUp & fVrDw
        IF (.NOT.implicitViscosity) THEN
@@ -403,7 +435,13 @@
 C-- No-slip and drag BCs appear as body forces in cell abutting topography 
         IF (no_slip_sides) THEN
 C-     No-slip BCs impose a drag at walls...
-         CALL MOM_U_SIDEDRAG(bi,bj,k,uFld,v4F,hFacZ,vF,myThid)
+         CALL MOM_U_SIDEDRAG(
+     I        bi,bj,k,
+     I        uFld, v4f, hFacZ,
+     I        viscAh_Z,viscA4_Z,
+     I        harmonic,biharmonic,useVariableViscosity,
+     O        vF,
+     I        myThid)
          DO j=jMin,jMax
           DO i=iMin,iMax
            gUdiss(i,j) = gUdiss(i,j) + vF(i,j)
@@ -542,11 +580,11 @@
 
 C     Laplacian and bi-harmonic terms, Zonal  Fluxes -> fZon
         CALL MOM_V_XVISCFLUX(bi,bj,k,vFld,v4f,hFacZ,fZon,
-     I    viscAh_Z,viscAh_D,viscA4_Z,viscA4_D,myThid)
+     I    viscAh_Z,viscA4_Z,myThid)
 
 C     Laplacian and bi-harmonic termis, Merid Fluxes -> fMer
         CALL MOM_V_YVISCFLUX(bi,bj,k,vFld,v4f,fMer,
-     I    viscAh_Z,viscAh_D,viscA4_Z,viscA4_D,myThid)
+     I    viscAh_D,viscA4_D,myThid)
 
 C     Eddy component of vertical flux (interior component only) -> fVrUp & fVrDw
        IF (.NOT.implicitViscosity) THEN
@@ -584,7 +622,13 @@
 C-- No-slip and drag BCs appear as body forces in cell abutting topography 
       IF (no_slip_sides) THEN
 C-     No-slip BCs impose a drag at walls...
-         CALL MOM_V_SIDEDRAG(bi,bj,k,vFld,v4F,hFacZ,vF,myThid)
+         CALL MOM_V_SIDEDRAG(
+     I        bi,bj,k,
+     I        vFld, v4f, hFacZ,
+     I        viscAh_Z,viscA4_Z,
+     I        harmonic,biharmonic,useVariableViscosity,
+     O        vF,
+     I        myThid)
          DO j=jMin,jMax
           DO i=iMin,iMax
            gvDiss(i,j) = gvDiss(i,j) + vF(i,j)
@@ -688,8 +732,7 @@
 
 #ifdef ALLOW_DIAGNOSTICS
       IF ( useDiagnostics ) THEN
-       IF (bottomDragTerms)
-     &  CALL DIAGNOSTICS_FILL(KE,    'momKE   ',k,1,2,bi,bj,myThid)
+        CALL DIAGNOSTICS_FILL(KE,    'momKE   ',k,1,2,bi,bj,myThid)
         CALL DIAGNOSTICS_FILL(gU(1-Olx,1-Oly,k,bi,bj),
      &                               'Um_Advec',k,1,2,bi,bj,myThid)
         CALL DIAGNOSTICS_FILL(gV(1-Olx,1-Oly,k,bi,bj),