--- MITgcm/model/src/dynamics.F	2005/05/23 20:49:37	1.117
+++ MITgcm/model/src/dynamics.F	2005/07/30 23:39:48	1.122
@@ -1,4 +1,4 @@
-C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/model/src/dynamics.F,v 1.117 2005/05/23 20:49:37 molod Exp $
+C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/model/src/dynamics.F,v 1.122 2005/07/30 23:39:48 jmc Exp $
 C $Name:  $
 
 #include "PACKAGES_CONFIG.h"
@@ -30,6 +30,7 @@
 C     | with the free-surface evolution or the rigid-lid:
 C     |   U[n] = U* + dt x d/dx P
 C     |   V[n] = V* + dt x d/dy P
+C     |   W[n] = W* + dt x d/dz P  (NH mode)
 C     |
 C     | "Calculation of Gs"
 C     | ===================
@@ -90,6 +91,8 @@
 C     !CALLING SEQUENCE:
 C     DYNAMICS()
 C      |
+C      |-- CALC_EP_FORCING
+C      |
 C      |-- CALC_GRAD_PHI_SURF
 C      |
 C      |-- CALC_VISCOSITY
@@ -104,11 +107,17 @@
 C      |
 C      |-- OBCS_APPLY_UV 
 C      |
+C      |-- MOM_U_IMPLICIT_R      
+C      |-- MOM_V_IMPLICIT_R      
+C      |
 C      |-- IMPLDIFF      
 C      |
 C      |-- OBCS_APPLY_UV
 C      |
-C      |-- CALL DEBUG_STATS_RL
+C      |-- CALC_GW
+C      |
+C      |-- DIAGNOSTICS_FILL
+C      |-- DEBUG_STATS_RL
 
 C     !INPUT/OUTPUT PARAMETERS:
 C     == Routine arguments ==
@@ -161,9 +170,7 @@
       INTEGER k, km1, kp1, kup, kDown
 
 #ifdef ALLOW_DIAGNOSTICS
-      _RL tmpFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
-      LOGICAL  DIAGNOSTICS_IS_ON
-      EXTERNAL DIAGNOSTICS_IS_ON
+      _RL tmpFac
 #endif /* ALLOW_DIAGNOSTICS */
 
  
@@ -264,8 +271,8 @@
 cph(
 c--   need some re-initialisation here to break dependencies
 cph)
-           gu(i,j,k,bi,bj) = 0. _d 0
-           gv(i,j,k,bi,bj) = 0. _d 0
+           gU(i,j,k,bi,bj) = 0. _d 0
+           gV(i,j,k,bi,bj) = 0. _d 0
 #endif
           ENDDO
          ENDDO
@@ -372,14 +379,15 @@
 #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,
+     I         KappaRU, KappaRV,
      U         fVerU, fVerV,
+     O         guDissip, gvDissip,
      I         myTime, myIter, myThid)
 #endif
 #ifdef ALLOW_MOM_VECINV
            IF (vectorInvariantMomentum) CALL MOM_VECINV(
      I         bi,bj,iMin,iMax,jMin,jMax,k,kup,kDown,
-     I         dPhiHydX,dPhiHydY,KappaRU,KappaRV,
+     I         KappaRU, KappaRV,
      U         fVerU, fVerV,
      O         guDissip, gvDissip,
      I         myTime, myIter, myThid)
@@ -476,6 +484,21 @@
 
 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
 
+#ifdef ALLOW_NONHYDROSTATIC
+C--   Step forward W field in N-H algorithm
+      IF ( momStepping .AND. nonHydrostatic ) THEN
+#ifdef ALLOW_DEBUG
+          IF ( debugLevel .GE. debLevB )
+     &     CALL DEBUG_CALL('CALC_GW',myThid)
+#endif
+         CALL TIMER_START('CALC_GW          [DYNAMICS]',myThid)
+         CALL CALC_GW( myTime, myIter, myThid )
+         CALL TIMER_STOP ('CALC_GW          [DYNAMICS]',myThid)
+      ENDIF
+#endif
+
+C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
+
 Cml(
 C     In order to compare the variance of phiHydLow of a p/z-coordinate
 C     run with etaH of a z/p-coordinate run the drift of phiHydLow
@@ -488,36 +511,14 @@
       IF ( usediagnostics ) THEN
 
        CALL DIAGNOSTICS_FILL(totPhihyd,'PHIHYD  ',0,Nr,0,1,1,myThid)
+       CALL DIAGNOSTICS_FILL(phiHydLow,'PHIBOT  ',0, 1,0,1,1,myThid)
 
-       IF ( DIAGNOSTICS_IS_ON('PHIHYDSQ',myThid) ) THEN
-        DO bj = myByLo(myThid), myByHi(myThid)
-        DO bi = myBxLo(myThid), myBxHi(myThid)
-         DO k = 1,Nr
-          DO j = 1,sNy
-           DO i = 1,sNx
-             tmpFld(i,j) = totPhihyd(i,j,k,bi,bj)*totPhihyd(i,j,k,bi,bj)
-           ENDDO
-          ENDDO
-          CALL DIAGNOSTICS_FILL(tmpFld,'PHIHYDSQ',k,1,2,bi,bj,myThid)
-         ENDDO
-        ENDDO
-        ENDDO
-       ENDIF
+       tmpFac = 1. _d 0
+       CALL DIAGNOSTICS_SCALE_FILL(totPhihyd,tmpFac,2,
+     &                                 'PHIHYDSQ',0,Nr,0,1,1,myThid)
 
-       CALL DIAGNOSTICS_FILL(phiHydLow,'PHIBOT  ',0,1,0,1,1,myThid)
-
-       IF ( DIAGNOSTICS_IS_ON('PHIBOTSQ',myThid) ) THEN
-        DO bj = myByLo(myThid), myByHi(myThid)
-         DO bi = myBxLo(myThid), myBxHi(myThid)
-          DO j = 1,sNy
-           DO i = 1,sNx
-             tmpFld(i,j) = phiHydLow(i,j,bi,bj)*phiHydLow(i,j,bi,bj)
-           ENDDO
-          ENDDO
-          CALL DIAGNOSTICS_FILL(tmpFld,'PHIBOTSQ',0,1,2,bi,bj,myThid)
-         ENDDO
-        ENDDO
-       ENDIF
+       CALL DIAGNOSTICS_SCALE_FILL(phiHydLow,tmpFac,2,
+     &                                 'PHIBOTSQ',0, 1,0,1,1,myThid)
 
       ENDIF
 #endif /* ALLOW_DIAGNOSTICS */