--- MITgcm/model/src/dynamics.F	2003/07/08 15:00:26	1.98
+++ MITgcm/model/src/dynamics.F	2004/07/06 01:05:53	1.107
@@ -1,6 +1,7 @@
-C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/model/src/dynamics.F,v 1.98 2003/07/08 15:00:26 heimbach Exp $
+C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/model/src/dynamics.F,v 1.107 2004/07/06 01:05:53 jmc Exp $
 C $Name:  $
 
+#include "PACKAGES_CONFIG.h"
 #include "CPP_OPTIONS.h"
 
 CBOP
@@ -72,6 +73,9 @@
 #include "EEPARAMS.h"
 #include "PARAMS.h"
 #include "DYNVARS.h"
+#ifdef ALLOW_CD_CODE
+#include "CD_CODE_VARS.h"
+#endif
 #include "GRID.h"
 #ifdef ALLOW_PASSIVE_TRACER
 #include "TR1.h"
@@ -320,6 +324,13 @@
        ENDDO
 #endif
 
+#ifdef ALLOW_AUTODIFF_TAMC
+CADJ STORE KappaRU(:,:,:) 
+CADJ &                 = comlev1_bibj, key=idynkey, byte=isbyte
+CADJ STORE KappaRV(:,:,:) 
+CADJ &                 = comlev1_bibj, key=idynkey, byte=isbyte
+#endif /* ALLOW_AUTODIFF_TAMC */
+
 C--     Start of dynamics loop
         DO k=1,Nr
 
@@ -334,40 +345,39 @@
 
 #ifdef ALLOW_AUTODIFF_TAMC 
          kkey = (idynkey-1)*Nr + k
+c
 CADJ STORE totphihyd (:,:,k,bi,bj) 
 CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
+CADJ STORE gt (:,:,k,bi,bj) 
+CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
+CADJ STORE gs (:,:,k,bi,bj) 
+CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
+CADJ STORE theta (:,:,k,bi,bj) 
+CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
+CADJ STORE salt  (:,:,k,bi,bj) 
+CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
 #endif /* ALLOW_AUTODIFF_TAMC */
 
 C--      Integrate hydrostatic balance for phiHyd with BC of 
 C        phiHyd(z=0)=0
-C        distinguishe between Stagger and Non Stagger time stepping
-         IF (staggerTimeStep) THEN
-           CALL CALC_PHI_HYD(
-     I        bi,bj,iMin,iMax,jMin,jMax,k,
-     I        gT, gS,
-     U        phiHydF,
-     O        phiHydC, dPhiHydX, dPhiHydY,
-     I        myTime, myIter, myThid )
-         ELSE
-           CALL CALC_PHI_HYD(
+         CALL CALC_PHI_HYD(
      I        bi,bj,iMin,iMax,jMin,jMax,k,
      I        theta, salt,
      U        phiHydF,
      O        phiHydC, dPhiHydX, dPhiHydY,
      I        myTime, myIter, myThid )
-         ENDIF
 
 C--      Calculate accelerations in the momentum equations (gU, gV, ...)
 C        and step forward storing the result in gU, gV, etc...
          IF ( momStepping ) THEN
-#ifndef DISABLE_MOM_FLUXFORM
+#ifdef ALLOW_MOM_FLUXFORM
            IF (.NOT. vectorInvariantMomentum) CALL MOM_FLUXFORM(
      I         bi,bj,iMin,iMax,jMin,jMax,k,kup,kDown,
      I         dPhiHydX,dPhiHydY,KappaRU,KappaRV,
      U         fVerU, fVerV,
      I         myTime, myIter, myThid)
 #endif
-#ifndef DISABLE_MOM_VECINV
+#ifdef ALLOW_MOM_VECINV
            IF (vectorInvariantMomentum) CALL MOM_VECINV(
      I         bi,bj,iMin,iMax,jMin,jMax,k,kup,kDown,
      I         dPhiHydX,dPhiHydY,KappaRU,KappaRV,
@@ -392,9 +402,19 @@
 C--     end of dynamics k loop (1:Nr)
         ENDDO
 
-C--     Implicit viscosity
-        IF (implicitViscosity.AND.momStepping) THEN
+C--     Implicit Vertical advection & viscosity
+#ifdef INCLUDE_IMPLVERTADV_CODE
+        IF ( momImplVertAdv ) THEN
+          CALL MOM_U_IMPLICIT_R( kappaRU, 
+     I                           bi, bj, myTime, myIter, myThid )
+          CALL MOM_V_IMPLICIT_R( kappaRV,
+     I                           bi, bj, myTime, myIter, myThid )
+        ELSEIF ( implicitViscosity ) THEN
+#else /* INCLUDE_IMPLVERTADV_CODE */
+        IF     ( implicitViscosity ) THEN
+#endif /* INCLUDE_IMPLVERTADV_CODE */
 #ifdef    ALLOW_AUTODIFF_TAMC
+CADJ STORE KappaRU(:,:,:) = comlev1_bibj , key=idynkey, byte=isbyte
 CADJ STORE gU(:,:,:,bi,bj) = comlev1_bibj , key=idynkey, byte=isbyte
 #endif    /* ALLOW_AUTODIFF_TAMC */
           CALL IMPLDIFF(
@@ -403,6 +423,7 @@
      U         gU,
      I         myThid )
 #ifdef    ALLOW_AUTODIFF_TAMC
+CADJ STORE KappaRV(:,:,:) = comlev1_bibj , key=idynkey, byte=isbyte
 CADJ STORE gV(:,:,:,bi,bj) = comlev1_bibj , key=idynkey, byte=isbyte
 #endif    /* ALLOW_AUTODIFF_TAMC */
           CALL IMPLDIFF(
@@ -410,17 +431,19 @@
      I         deltaTmom, KappaRV,recip_HFacS,
      U         gV,
      I         myThid )
+        ENDIF
 
 #ifdef   ALLOW_OBCS
 C--      Apply open boundary conditions
-         IF (useOBCS) THEN
+        IF ( useOBCS .AND.(implicitViscosity.OR.momImplVertAdv) ) THEN
            DO K=1,Nr
              CALL OBCS_APPLY_UV( bi, bj, k, gU, gV, myThid )
            ENDDO
-         END IF
+        ENDIF
 #endif   /* ALLOW_OBCS */
 
-#ifdef    INCLUDE_CD_CODE
+#ifdef    ALLOW_CD_CODE
+        IF (implicitViscosity.AND.useCDscheme) THEN
 #ifdef    ALLOW_AUTODIFF_TAMC
 CADJ STORE vVelD(:,:,:,bi,bj) = comlev1_bibj , key=idynkey, byte=isbyte
 #endif    /* ALLOW_AUTODIFF_TAMC */
@@ -437,9 +460,9 @@
      I         deltaTmom, KappaRV,recip_HFacS,
      U         uVelD,
      I         myThid )
-#endif    /* INCLUDE_CD_CODE */
-C--     End If implicitViscosity.AND.momStepping
         ENDIF
+#endif    /* ALLOW_CD_CODE */
+C--     End implicit Vertical advection & viscosity
  
        ENDDO
       ENDDO
@@ -452,7 +475,7 @@
 C     &                'phiHydLow', myThid )
 Cml)
 
-#ifndef DISABLE_DEBUGMODE
+#ifdef ALLOW_DEBUG
       If ( debugLevel .GE. debLevB ) THEN
        CALL DEBUG_STATS_RL(1,EtaN,'EtaN (DYNAMICS)',myThid)
        CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (DYNAMICS)',myThid)