--- MITgcm/model/src/dynamics.F	1998/11/06 22:44:45	1.38
+++ MITgcm/model/src/dynamics.F	1999/05/18 18:01:12	1.42
@@ -1,4 +1,4 @@
-C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/model/src/dynamics.F,v 1.38 1998/11/06 22:44:45 cnh Exp $
+C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/model/src/dynamics.F,v 1.42 1999/05/18 18:01:12 adcroft Exp $
 
 #include "CPP_OPTIONS.h"
 
@@ -20,6 +20,7 @@
 C     | C*P* comments indicating place holders for which code is |
 C     |      presently being developed.                          |
 C     \==========================================================/
+      IMPLICIT NONE
 
 C     == Global variables ===
 #include "SIZE.h"
@@ -27,6 +28,10 @@
 #include "CG2D.h"
 #include "PARAMS.h"
 #include "DYNVARS.h"
+#include "GRID.h"
+#ifdef ALLOW_KPP
+#include "KPPMIX.h"
+#endif
 
 C     == Routine arguments ==
 C     myTime - Current time in simulation
@@ -117,6 +122,8 @@
       _RL KapGM   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
       _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
       _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
+      _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
+      _RL KappaRV (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
 
       INTEGER iMin, iMax
       INTEGER jMin, jMax
@@ -251,11 +258,13 @@
         CALL CORRECTION_STEP(
      I       bi,bj,iMin,iMax,jMin,jMax,K,
      I       etaSurfX,etaSurfY,myTime,myThid)
+        IF (openBoundaries) CALL APPLY_OBCS1( bi, bj, K, myThid )
         IF ( .NOT. BOTTOM_LAYER ) THEN
 C--      Update fields in layer below according to tendency terms
          CALL CORRECTION_STEP(
      I        bi,bj,iMin,iMax,jMin,jMax,K+1,
      I        etaSurfX,etaSurfY,myTime,myThid)
+         IF (openBoundaries) CALL APPLY_OBCS1( bi, bj, K+1, myThid )
         ENDIF
 #endif
 C--     Density of 1st level (below W(1)) reference to level 1
@@ -266,7 +275,11 @@
      I     myThid )
 #endif
 
-        IF ( .NOT. BOTTOM_LAYER ) THEN
+        IF (       (.NOT. BOTTOM_LAYER)
+#ifdef ALLOW_KPP
+     &       .AND. (.NOT.usingKPPmixing) ! CONVECT not needed with KPP mixing
+#endif
+     &     ) THEN
 C--      Check static stability with layer below
 C--      and mix as needed.
 #ifdef  INCLUDE_FIND_RHO_CALL
@@ -308,6 +321,7 @@
           CALL CORRECTION_STEP(
      I         bi,bj,iMin,iMax,jMin,jMax,K+1,
      I         etaSurfX,etaSurfY,myTime,myThid)
+          IF (openBoundaries) CALL APPLY_OBCS1( bi, bj, K+1, myThid )
          ENDIF
 #endif
 C--      Density of K level (below W(K)) reference to K level
@@ -317,7 +331,11 @@
      O      rhoK,
      I      myThid )
 #endif
-         IF ( .NOT. BOTTOM_LAYER ) THEN
+         IF (       (.NOT. BOTTOM_LAYER)
+#ifdef ALLOW_KPP
+     &       .AND. (.NOT.usingKPPmixing) ! CONVECT not needed with KPP mixing
+#endif
+     &      ) THEN
 C--       Check static stability with layer below and mix as needed.
 C--       Density of K+1 level (below W(K+1)) reference to K level.
 #ifdef  INCLUDE_FIND_RHO_CALL
@@ -374,6 +392,18 @@
          ENDDO
         ENDDO ! K
 
+#ifdef ALLOW_KPP
+C--     Compute KPP mixing coefficients
+        IF (usingKPPmixing) THEN
+         CALL TIMER_START('KVMIX (FIND KPP COEFFICIENTS) [DYNAMICS]'
+     I          , myThid)
+         CALL KVMIX(
+     I               bi, bj, myTime, myThid )
+         CALL TIMER_STOP ('KVMIX (FIND KPP COEFFICIENTS) [DYNAMICS]'
+     I        , myThid)
+        ENDIF
+#endif
+
         DO K = Nr, 1, -1
 
          kM1  =max(1,k-1)   ! Points to level above k (=k-1)
@@ -394,7 +424,7 @@
          CALL CALC_DIFFUSIVITY(
      I        bi,bj,iMin,iMax,jMin,jMax,K,
      I        maskC,maskUp,KapGM,K33,
-     O        KappaRT,KappaRS,
+     O        KappaRT,KappaRS,KappaRU,KappaRV,
      I        myThid)
 #endif
 C--      Calculate accelerations in the momentum equations
@@ -402,7 +432,7 @@
           CALL CALC_MOM_RHS(
      I         bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,
      I         xA,yA,uTrans,vTrans,rTrans,rVel,maskC,
-     I         phiHyd,
+     I         phiHyd,KappaRU,KappaRV,
      U         aTerm,xTerm,cTerm,mTerm,pTerm,
      U         fZon, fMer, fVerU, fVerV,
      I         myTime, myThid)
@@ -428,6 +458,11 @@
          CALL TIMESTEP(
      I       bi,bj,iMin,iMax,jMin,jMax,K,
      I       myThid)
+C--      Apply open boundary conditions
+         IF (openBoundaries) CALL APPLY_OBCS2( bi, bj, K, myThid )
+C--      Freeze water
+         IF (allowFreezing)
+     &   CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, K, myThid )
 C--      Diagnose barotropic divergence of predicted fields
          CALL CALC_DIV_GHAT(
      I       bi,bj,iMin,iMax,jMin,jMax,K,
@@ -448,10 +483,41 @@
 
 C--     Implicit diffusion
         IF (implicitDiffusion) THEN
-         CALL IMPLDIFF( bi, bj, iMin, iMax, jMin, jMax,
-     I                  KappaRT,KappaRS,
-     I                  myThid )
-        ENDIF
+         IF (tempStepping) CALL IMPLDIFF(
+     I         bi, bj, iMin, iMax, jMin, jMax,
+     I         deltaTtracer, KappaRT,recip_HFacC,
+     U         gTNm1,
+     I         myThid )
+         IF (saltStepping) CALL IMPLDIFF(
+     I         bi, bj, iMin, iMax, jMin, jMax,
+     I         deltaTtracer, KappaRS,recip_HFacC,
+     U         gSNm1,
+     I         myThid )
+         IF (momStepping) THEN
+          CALL IMPLDIFF(
+     I         bi, bj, iMin, iMax, jMin, jMax,
+     I         deltaTmom, KappaRU,recip_HFacW,
+     U         gUNm1,
+     I         myThid )
+          CALL IMPLDIFF(
+     I         bi, bj, iMin, iMax, jMin, jMax,
+     I         deltaTmom, KappaRV,recip_HFacS,
+     U         gVNm1,
+     I         myThid )
+#ifdef INCLUDE_CD_CODE
+          CALL IMPLDIFF(
+     I         bi, bj, iMin, iMax, jMin, jMax,
+     I         deltaTmom, KappaRU,recip_HFacW,
+     U         vVelD,
+     I         myThid )
+          CALL IMPLDIFF(
+     I         bi, bj, iMin, iMax, jMin, jMax,
+     I         deltaTmom, KappaRV,recip_HFacS,
+     U         uVelD,
+     I         myThid )
+#endif
+         ENDIF ! momStepping
+        ENDIF ! implicitDiffusion
  
        ENDDO
       ENDDO