--- MITgcm/pkg/mom_fluxform/mom_fluxform.F	2003/10/11 16:37:55	1.15
+++ MITgcm/pkg/mom_fluxform/mom_fluxform.F	2004/10/21 21:14:02	1.20
@@ -1,4 +1,4 @@
-C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/mom_fluxform/mom_fluxform.F,v 1.15 2003/10/11 16:37:55 jmc Exp $
+C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/mom_fluxform/mom_fluxform.F,v 1.20 2004/10/21 21:14:02 afe Exp $
 C $Name:  $
 
 CBOI
@@ -122,6 +122,15 @@
       _RL  vFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
       _RL  rTransU(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
       _RL  rTransV(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
+      _RL KE(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
+      _RL viscAhD(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
+      _RL viscAhZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
+      _RL viscA4D(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
+      _RL viscA4Z(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
+      _RL vort3(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
+      _RL hDiv(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
+      _RL strain(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
+      _RL tension(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
 C     I,J,K - Loop counters
 C     rVelMaskOverride - Factor for imposing special surface boundary conditions
 C                        ( set according to free-surface condition ).
@@ -155,7 +164,6 @@
       INTEGER km1,kp1
       _RL wVelBottomOverride
       LOGICAL bottomDragTerms
-      _RL KE(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
 CEOP
 
       km1=MAX(1,k-1)
@@ -179,14 +187,8 @@
         fMer(i,j) = 0.
         rTransU(i,j) = 0.
         rTransV(i,j) = 0.
-#ifdef ALLOW_AUTODIFF_TAMC
-C- jmc: this is wrong, but at least with #ifdef/endif TAMC, it does not break 
-C       the forward code ; (same thing in mom_vectinv)
-        fVerU(i,j,1) = 0. _d 0
-        fVerU(i,j,2) = 0. _d 0
-        fVerV(i,j,1) = 0. _d 0
-        fVerV(i,j,2) = 0. _d 0
-#endif
+        strain(i,j) = 0.
+        tension(i,j) = 0.
        ENDDO
       ENDDO
 
@@ -261,7 +263,16 @@
        ENDDO
       ENDDO
 
-      CALL MOM_CALC_KE(bi,bj,k,uFld,vFld,KE,myThid)
+      CALL MOM_CALC_KE(bi,bj,k,3,uFld,vFld,KE,myThid)
+
+c     IF (viscAstrain.NE.0. .OR. viscAtension.NE.0.) THEN
+         CALL MOM_CALC_TENSION(bi,bj,k,uFld,vFld,
+     O                         tension,
+     I                         myThid)
+         CALL MOM_CALC_STRAIN(bi,bj,k,uFld,vFld,hFacZ,
+     O                        strain,
+     I                        myThid)
+c     ENDIF
 
 C---  First call (k=1): compute vertical adv. flux fVerU(kUp) & fVerV(kUp)
       IF (momAdvection.AND.k.EQ.1) THEN
@@ -297,6 +308,10 @@
      I                       myTime, myIter, myThid)
       ENDIF
 
+c     IF (momViscosity) THEN
+c    &  CALL MOM_CALC_VISCOSITY(bi,bj,k,
+c    I                         uFld,vFld,
+c    O                         viscAhD,viscAhZ,myThid)
 
 C---- Zonal momentum equation starts here
 
@@ -440,8 +455,17 @@
          gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mTFacU*mT(i,j)
         ENDDO
        ENDDO
+                                                                                
+      ENDIF
+      IF (usingCylindricalGrid) THEN
+         CALL MOM_U_METRIC_CYLINDER(bi,bj,k,uFld,vFld,mT,myThid)
+         DO j=jMin,jMax
+          DO i=iMin,iMax
+             gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)+mTFacU*mT(i,j)
+          ENDDO
+       ENDDO
+                                                                                
       ENDIF
-
 C--   Set du/dt on boundaries to zero
       DO j=jMin,jMax
        DO i=iMin,iMax
@@ -593,6 +617,14 @@
         ENDDO
        ENDDO
       ENDIF
+      IF (usingCylindricalGrid) THEN
+         CALL MOM_V_METRIC_CYLINDER(bi,bj,k,uFld,vFld,mT,myThid)
+         DO j=jMin,jMax
+            DO i=iMin,iMax
+               gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)+mTFacV*mT(i,j)
+            ENDDO
+         ENDDO
+      ENDIF
 
 C--   Set dv/dt on boundaries to zero
       DO j=jMin,jMax