--- MITgcm/model/src/dynamics.F	2001/02/07 21:48:02	1.61
+++ MITgcm/model/src/dynamics.F	2001/03/06 16:59:44	1.64
@@ -1,4 +1,4 @@
-C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/model/src/dynamics.F,v 1.61 2001/02/07 21:48:02 jmc Exp $
+C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/model/src/dynamics.F,v 1.64 2001/03/06 16:59:44 jmc Exp $
 C $Name:  $
 
 #include "CPP_OPTIONS.h"
@@ -26,7 +26,6 @@
 C     == Global variables ===
 #include "SIZE.h"
 #include "EEPARAMS.h"
-#include "CG2D.h"
 #include "PARAMS.h"
 #include "DYNVARS.h"
 #include "GRID.h"
@@ -40,6 +39,10 @@
 # include "KPP.h"
 #endif
 
+#ifdef ALLOW_TIMEAVE
+#include "TIMEAVE_STATV.h"
+#endif
+
 C     == Routine arguments ==
 C     myTime - Current time in simulation
 C     myIter - Current iteration number in simulation
@@ -68,8 +71,8 @@
 C                      In p coords phiHydiHyd is the geopotential 
 C                      surface height 
 C                      anomaly.
-C     etaSurfX,      - Holds surface elevation gradient in X and Y.
-C     etaSurfY
+C     phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean)
+C     phiSurfY             or geopotentiel (atmos) in X and Y direction
 C     KappaRT,       - Total diffusion in vertical for T and S.
 C     KappaRS          (background + spatially varying, isopycnal term).
 C     iMin, iMax     - Ranges and sub-block indices on which calculations
@@ -92,6 +95,8 @@
       _RL phiHyd  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
       _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
       _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
+      _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
+      _RL phiSurfY(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)
@@ -100,10 +105,8 @@
       _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
       _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
 
-C This is currently also used by IVDC and Diagnostics
-C #ifdef INCLUDE_CONVECT_CALL
+C This is currently used by IVDC and Diagnostics
       _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
-C #endif
 
       INTEGER iMin, iMax
       INTEGER jMin, jMax
@@ -111,6 +114,12 @@
       INTEGER i, j
       INTEGER k, km1, kup, kDown
 
+Cjmc : add for phiHyd output <- but not working if multi tile per CPU
+c     CHARACTER*(MAX_LEN_MBUF) suff 
+c     LOGICAL  DIFFERENT_MULTIPLE
+c     EXTERNAL DIFFERENT_MULTIPLE
+Cjmc(end)
+ 
 #ifdef ALLOW_AUTODIFF_TAMC
       INTEGER    isbyte
       PARAMETER( isbyte = 4 )
@@ -192,6 +201,8 @@
         rhoKM1 (i,j) = 0. _d 0
         rhok   (i,j) = 0. _d 0
         maskC  (i,j) = 0. _d 0
+        phiSurfX(i,j) = 0. _d 0
+        phiSurfY(i,j) = 0. _d 0
        ENDDO
       ENDDO
 
@@ -248,9 +259,8 @@
         DO k=1,Nr
          DO j=1-OLy,sNy+OLy
           DO i=1-OLx,sNx+OLx
-#ifdef INCLUDE_CONVECT_CALL
+C This is currently also used by IVDC and Diagnostics
            ConvectCount(i,j,k) = 0.
-#endif
            KappaRT(i,j,k) = 0. _d 0
            KappaRS(i,j,k) = 0. _d 0
           ENDDO
@@ -316,7 +326,7 @@
      I        rhoKm1, rhoK,
      U        ConvectCount, KappaRT, KappaRS,
      I        myTime, myIter, myThid)
-          END IF
+          ENDIF
 
 C--     end of diagnostic k loop (Nr:1)
         ENDDO
@@ -482,7 +492,7 @@
 C--     Implicit diffusion
         IF (implicitDiffusion) THEN
 
-          IF (tempStepping) THEN
+         IF (tempStepping) THEN
 #ifdef ALLOW_AUTODIFF_TAMC
             idkey = iikey + 1
 #endif /* ALLOW_AUTODIFF_TAMC */
@@ -516,7 +526,24 @@
 C--     End If implicitDiffusion
         ENDIF
 
-
+C--     Start computation of dynamics
+        iMin = 1-OLx+2
+        iMax = sNx+OLx-1
+        jMin = 1-OLy+2
+        jMax = sNy+OLy-1
+
+C--     Explicit part of the Surface Pressure Gradient (add in TIMESTEP)
+C       (note: this loop will be replaced by CALL CALC_GRAD_ETA)
+        IF (implicSurfPress.NE.1.) THEN
+          DO j=jMin,jMax
+            DO i=iMin,iMax 
+              phiSurfX(i,j) = _recip_dxC(i,j,bi,bj)*gBaro
+     &           *(etaN(i,j,bi,bj)-etaN(i-1,j,bi,bj)) 
+              phiSurfY(i,j) = _recip_dyC(i,j,bi,bj)*gBaro 
+     &           *(etaN(i,j,bi,bj)-etaN(i,j-1,bi,bj)) 
+            ENDDO
+          ENDDO
+        ENDIF
 
 C--     Start of dynamics loop
         DO k=1,Nr
@@ -529,11 +556,6 @@
           kup  = 1+MOD(k+1,2)
           kDown= 1+MOD(k,2)
 
-          iMin = 1-OLx+2
-          iMax = sNx+OLx-1
-          jMin = 1-OLy+2
-          jMax = sNy+OLy-1
-
 C--      Integrate hydrostatic balance for phiHyd with BC of 
 C        phiHyd(z=0)=0
 C        distinguishe between Stagger and Non Stagger time stepping
@@ -560,7 +582,8 @@
      U         fVerU, fVerV,
      I         myTime, myThid)
            CALL TIMESTEP(
-     I         bi,bj,iMin,iMax,jMin,jMax,k,phiHyd,
+     I         bi,bj,iMin,iMax,jMin,jMax,k,
+     I         phiHyd, phiSurfX, phiSurfY,
      I         myIter, myThid)
 
 #ifdef   ALLOW_OBCS
@@ -638,6 +661,25 @@
 C--     End If implicitViscosity.AND.momStepping
         ENDIF
  
+Cjmc : add for phiHyd output <- but not working if multi tile per CPU
+c       IF ( DIFFERENT_MULTIPLE(dumpFreq,myTime+deltaTClock,myTime) 
+c    &  .AND. buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN
+c         WRITE(suff,'(I10.10)') myIter+1
+c         CALL WRITE_FLD_XYZ_RL('PH.',suff,phiHyd,myIter+1,myThid) 
+c       ENDIF
+Cjmc(end)
+
+#ifdef ALLOW_TIMEAVE
+        IF (taveFreq.GT.0.) THEN 
+          CALL TIMEAVE_CUMULATE(phiHydtave, phiHyd, Nr,
+     I                              deltaTclock, bi, bj, myThid)
+          IF (ivdc_kappa.NE.0.) THEN
+            CALL TIMEAVE_CUMULATE(ConvectCountTave, ConvectCount, Nr,
+     I                              deltaTclock, bi, bj, myThid)
+          ENDIF
+        ENDIF
+#endif /* ALLOW_TIMEAVE */
+
        ENDDO
       ENDDO