/[MITgcm]/MITgcm/model/src/dynamics.F
ViewVC logotype

Diff of /MITgcm/model/src/dynamics.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.74 by heimbach, Mon Jul 30 20:37:45 2001 UTC revision 1.90 by mlosch, Wed Sep 18 16:38:01 2002 UTC
# Line 3  C $Name$ Line 3  C $Name$
3    
4  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
5    
6    CBOP
7    C     !ROUTINE: DYNAMICS
8    C     !INTERFACE:
9        SUBROUTINE DYNAMICS(myTime, myIter, myThid)        SUBROUTINE DYNAMICS(myTime, myIter, myThid)
10  C     /==========================================================\  C     !DESCRIPTION: \bv
11  C     | SUBROUTINE DYNAMICS                                      |  C     *==========================================================*
12  C     | o Controlling routine for the explicit part of the model |  C     | SUBROUTINE DYNAMICS                                      
13  C     |   dynamics.                                              |  C     | o Controlling routine for the explicit part of the model  
14  C     |==========================================================|  C     |   dynamics.                                              
15  C     | This routine evaluates the "dynamics" terms for each     |  C     *==========================================================*
16  C     | block of ocean in turn. Because the blocks of ocean have |  C     | This routine evaluates the "dynamics" terms for each      
17  C     | overlap regions they are independent of one another.     |  C     | block of ocean in turn. Because the blocks of ocean have  
18  C     | If terms involving lateral integrals are needed in this  |  C     | overlap regions they are independent of one another.      
19  C     | routine care will be needed. Similarly finite-difference |  C     | If terms involving lateral integrals are needed in this  
20  C     | operations with stencils wider than the overlap region   |  C     | routine care will be needed. Similarly finite-difference  
21  C     | require special consideration.                           |  C     | operations with stencils wider than the overlap region    
22  C     | Notes                                                    |  C     | require special consideration.                            
23  C     | =====                                                    |  C     | The algorithm...
24  C     | C*P* comments indicating place holders for which code is |  C     |
25  C     |      presently being developed.                          |  C     | "Correction Step"
26  C     \==========================================================/  C     | =================
27    C     | Here we update the horizontal velocities with the surface
28    C     | pressure such that the resulting flow is either consistent
29    C     | with the free-surface evolution or the rigid-lid:
30    C     |   U[n] = U* + dt x d/dx P
31    C     |   V[n] = V* + dt x d/dy P
32    C     |
33    C     | "Calculation of Gs"
34    C     | ===================
35    C     | This is where all the accelerations and tendencies (ie.
36    C     | physics, parameterizations etc...) are calculated
37    C     |   rho = rho ( theta[n], salt[n] )
38    C     |   b   = b(rho, theta)
39    C     |   K31 = K31 ( rho )
40    C     |   Gu[n] = Gu( u[n], v[n], wVel, b, ... )
41    C     |   Gv[n] = Gv( u[n], v[n], wVel, b, ... )
42    C     |   Gt[n] = Gt( theta[n], u[n], v[n], wVel, K31, ... )
43    C     |   Gs[n] = Gs( salt[n], u[n], v[n], wVel, K31, ... )
44    C     |
45    C     | "Time-stepping" or "Prediction"
46    C     | ================================
47    C     | The models variables are stepped forward with the appropriate
48    C     | time-stepping scheme (currently we use Adams-Bashforth II)
49    C     | - For momentum, the result is always *only* a "prediction"
50    C     | in that the flow may be divergent and will be "corrected"
51    C     | later with a surface pressure gradient.
52    C     | - Normally for tracers the result is the new field at time
53    C     | level [n+1} *BUT* in the case of implicit diffusion the result
54    C     | is also *only* a prediction.
55    C     | - We denote "predictors" with an asterisk (*).
56    C     |   U* = U[n] + dt x ( 3/2 Gu[n] - 1/2 Gu[n-1] )
57    C     |   V* = V[n] + dt x ( 3/2 Gv[n] - 1/2 Gv[n-1] )
58    C     |   theta[n+1] = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
59    C     |   salt[n+1] = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
60    C     | With implicit diffusion:
61    C     |   theta* = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
62    C     |   salt* = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
63    C     |   (1 + dt * K * d_zz) theta[n] = theta*
64    C     |   (1 + dt * K * d_zz) salt[n] = salt*
65    C     |
66    C     *==========================================================*
67    C     \ev
68    C     !USES:
69        IMPLICIT NONE        IMPLICIT NONE
   
70  C     == Global variables ===  C     == Global variables ===
71  #include "SIZE.h"  #include "SIZE.h"
72  #include "EEPARAMS.h"  #include "EEPARAMS.h"
# Line 32  C     == Global variables === Line 76  C     == Global variables ===
76  #ifdef ALLOW_PASSIVE_TRACER  #ifdef ALLOW_PASSIVE_TRACER
77  #include "TR1.h"  #include "TR1.h"
78  #endif  #endif
   
79  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
80  # include "tamc.h"  # include "tamc.h"
81  # include "tamc_keys.h"  # include "tamc_keys.h"
# Line 40  C     == Global variables === Line 83  C     == Global variables ===
83  # ifdef ALLOW_KPP  # ifdef ALLOW_KPP
84  #  include "KPP.h"  #  include "KPP.h"
85  # endif  # endif
 # ifdef ALLOW_GMREDI  
 #  include "GMREDI.h"  
 # endif  
86  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
   
87  #ifdef ALLOW_TIMEAVE  #ifdef ALLOW_TIMEAVE
88  #include "TIMEAVE_STATV.h"  #include "TIMEAVE_STATV.h"
89  #endif  #endif
90    
91    C     !CALLING SEQUENCE:
92    C     DYNAMICS()
93    C      |
94    C      |-- CALC_GRAD_PHI_SURF
95    C      |
96    C      |-- CALC_VISCOSITY
97    C      |
98    C      |-- CALC_PHI_HYD  
99    C      |
100    C      |-- STORE_PRESSURE
101    C      |
102    C      |-- MOM_FLUXFORM  
103    C      |
104    C      |-- MOM_VECINV    
105    C      |
106    C      |-- TIMESTEP      
107    C      |
108    C      |-- OBCS_APPLY_UV
109    C      |
110    C      |-- IMPLDIFF      
111    C      |
112    C      |-- OBCS_APPLY_UV
113    C      |
114    C      |-- CALL TIMEAVE_CUMUL_1T
115    C      |-- CALL DEBUG_STATS_RL
116    
117    C     !INPUT/OUTPUT PARAMETERS:
118  C     == Routine arguments ==  C     == Routine arguments ==
119  C     myTime - Current time in simulation  C     myTime - Current time in simulation
120  C     myIter - Current iteration number in simulation  C     myIter - Current iteration number in simulation
# Line 57  C     myThid - Thread number for this in Line 123  C     myThid - Thread number for this in
123        INTEGER myIter        INTEGER myIter
124        INTEGER myThid        INTEGER myThid
125    
126    C     !LOCAL VARIABLES:
127  C     == Local variables  C     == Local variables
 C     xA, yA                 - Per block temporaries holding face areas  
 C     uTrans, vTrans, rTrans - Per block temporaries holding flow  
 C                              transport  
 C                              o uTrans: Zonal transport  
 C                              o vTrans: Meridional transport  
 C                              o rTrans: Vertical transport  
 C     maskUp                   o maskUp: land/water mask for W points  
128  C     fVer[STUV]               o fVer: Vertical flux term - note fVer  C     fVer[STUV]               o fVer: Vertical flux term - note fVer
129  C                                      is "pipelined" in the vertical  C                                      is "pipelined" in the vertical
130  C                                      so we need an fVer for each  C                                      so we need an fVer for each
# Line 77  C                      In p coords phiHy Line 137  C                      In p coords phiHy
137  C                      surface height anomaly.  C                      surface height anomaly.
138  C     phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean)  C     phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean)
139  C     phiSurfY             or geopotentiel (atmos) in X and Y direction  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).  
140  C     iMin, iMax     - Ranges and sub-block indices on which calculations  C     iMin, iMax     - Ranges and sub-block indices on which calculations
141  C     jMin, jMax       are applied.  C     jMin, jMax       are applied.
142  C     bi, bj  C     bi, bj
143  C     k, kup,        - Index for layer above and below. kup and kDown  C     k, kup,        - Index for layer above and below. kup and kDown
144  C     kDown, km1       are switched with layer to be the appropriate  C     kDown, km1       are switched with layer to be the appropriate
145  C                      index into fVerTerm.  C                      index into fVerTerm.
 C     tauAB - Adams-Bashforth timestepping weight: 0=forward ; 1/2=Adams-Bashf.  
       _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)  
       _RL fVerS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)  
       _RL fVerTr1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)  
146        _RL fVerU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
147        _RL fVerV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
148        _RL phiHyd  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL phiHyd  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
# Line 102  C     tauAB - Adams-Bashforth timesteppi Line 150  C     tauAB - Adams-Bashforth timesteppi
150        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
151        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
152        _RL phiSurfY(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)  
153        _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
154        _RL KappaRV (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRV (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
       _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
       _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
       _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
       _RL tauAB  
   
 C This is currently used by IVDC and Diagnostics  
       _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
155    
156        INTEGER iMin, iMax        INTEGER iMin, iMax
157        INTEGER jMin, jMax        INTEGER jMin, jMax
158        INTEGER bi, bj        INTEGER bi, bj
159        INTEGER i, j        INTEGER i, j
160        INTEGER k, km1, kup, kDown        INTEGER k, km1, kp1, kup, kDown
161    
162  Cjmc : add for phiHyd output <- but not working if multi tile per CPU  Cjmc : add for phiHyd output <- but not working if multi tile per CPU
163  c     CHARACTER*(MAX_LEN_MBUF) suff  c      CHARACTER*(MAX_LEN_MBUF) suff
164  c     LOGICAL  DIFFERENT_MULTIPLE  c      LOGICAL  DIFFERENT_MULTIPLE
165  c     EXTERNAL DIFFERENT_MULTIPLE  c      EXTERNAL DIFFERENT_MULTIPLE
166  Cjmc(end)  Cjmc(end)
167    
168  C---    The algorithm...  C---    The algorithm...
# Line 169  C         salt* = salt[n] + dt x ( 3/2 G Line 208  C         salt* = salt[n] + dt x ( 3/2 G
208  C         (1 + dt * K * d_zz) theta[n] = theta*  C         (1 + dt * K * d_zz) theta[n] = theta*
209  C         (1 + dt * K * d_zz) salt[n] = salt*  C         (1 + dt * K * d_zz) salt[n] = salt*
210  C---  C---
211    CEOP
 #ifdef ALLOW_AUTODIFF_TAMC  
 C--   dummy statement to end declaration part  
       ikey = 1  
 #endif /* ALLOW_AUTODIFF_TAMC */  
212    
213  C--   Set up work arrays with valid (i.e. not NaN) values  C--   Set up work arrays with valid (i.e. not NaN) values
214  C     These inital values do not alter the numerical results. They  C     These inital values do not alter the numerical results. They
# Line 182  C     point numbers. This prevents spuri Line 217  C     point numbers. This prevents spuri
217  C     uninitialised but inert locations.  C     uninitialised but inert locations.
218        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
219         DO i=1-OLx,sNx+OLx         DO i=1-OLx,sNx+OLx
         xA(i,j)      = 0. _d 0  
         yA(i,j)      = 0. _d 0  
         uTrans(i,j)  = 0. _d 0  
         vTrans(i,j)  = 0. _d 0  
         DO k=1,Nr  
          phiHyd(i,j,k)  = 0. _d 0  
          KappaRU(i,j,k) = 0. _d 0  
          KappaRV(i,j,k) = 0. _d 0  
          sigmaX(i,j,k) = 0. _d 0  
          sigmaY(i,j,k) = 0. _d 0  
          sigmaR(i,j,k) = 0. _d 0  
         ENDDO  
220          rhoKM1 (i,j) = 0. _d 0          rhoKM1 (i,j) = 0. _d 0
221          rhok   (i,j) = 0. _d 0          rhok   (i,j) = 0. _d 0
222          phiSurfX(i,j) = 0. _d 0          phiSurfX(i,j) = 0. _d 0
# Line 201  C     uninitialised but inert locations. Line 224  C     uninitialised but inert locations.
224         ENDDO         ENDDO
225        ENDDO        ENDDO
226    
227    C-- Call to routine for calculation of
228    C   Eliassen-Palm-flux-forced U-tendency,
229    C   if desired:
230    #ifdef INCLUDE_EP_FORCING_CODE
231          CALL CALC_EP_FORCING(myThid)
232    #endif
233    
234  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
235  C--   HPF directive to help TAMC  C--   HPF directive to help TAMC
# Line 211  CHPF$ INDEPENDENT Line 240  CHPF$ INDEPENDENT
240    
241  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
242  C--    HPF directive to help TAMC  C--    HPF directive to help TAMC
243  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS,fVerU,fVerV  CHPF$  INDEPENDENT, NEW (fVerU,fVerV
244  CHPF$&                  ,phiHyd,utrans,vtrans,xA,yA  CHPF$&                  ,phiHyd
245  CHPF$&                  ,KappaRT,KappaRS,KappaRU,KappaRV  CHPF$&                  ,KappaRU,KappaRV
246  CHPF$&                  )  CHPF$&                  )
247  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
248    
# Line 222  CHPF$&                  ) Line 251  CHPF$&                  )
251  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
252            act1 = bi - myBxLo(myThid)            act1 = bi - myBxLo(myThid)
253            max1 = myBxHi(myThid) - myBxLo(myThid) + 1            max1 = myBxHi(myThid) - myBxLo(myThid) + 1
   
254            act2 = bj - myByLo(myThid)            act2 = bj - myByLo(myThid)
255            max2 = myByHi(myThid) - myByLo(myThid) + 1            max2 = myByHi(myThid) - myByLo(myThid) + 1
   
256            act3 = myThid - 1            act3 = myThid - 1
257            max3 = nTx*nTy            max3 = nTx*nTy
   
258            act4 = ikey_dynamics - 1            act4 = ikey_dynamics - 1
   
259            ikey = (act1 + 1) + act2*max1            ikey = (act1 + 1) + act2*max1
260       &                      + act3*max1*max2       &                      + act3*max1*max2
261       &                      + act4*max1*max2*max3       &                      + act4*max1*max2*max3
# Line 239  CHPF$&                  ) Line 264  CHPF$&                  )
264  C--     Set up work arrays that need valid initial values  C--     Set up work arrays that need valid initial values
265          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
266           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
267            rTrans (i,j)   = 0. _d 0            DO k=1,Nr
268            fVerT  (i,j,1) = 0. _d 0             phiHyd(i,j,k)  = 0. _d 0
269            fVerT  (i,j,2) = 0. _d 0             KappaRU(i,j,k) = 0. _d 0
270            fVerS  (i,j,1) = 0. _d 0             KappaRV(i,j,k) = 0. _d 0
271            fVerS  (i,j,2) = 0. _d 0            ENDDO
           fVerTr1(i,j,1) = 0. _d 0  
           fVerTr1(i,j,2) = 0. _d 0  
272            fVerU  (i,j,1) = 0. _d 0            fVerU  (i,j,1) = 0. _d 0
273            fVerU  (i,j,2) = 0. _d 0            fVerU  (i,j,2) = 0. _d 0
274            fVerV  (i,j,1) = 0. _d 0            fVerV  (i,j,1) = 0. _d 0
# Line 253  C--     Set up work arrays that need val Line 276  C--     Set up work arrays that need val
276           ENDDO           ENDDO
277          ENDDO          ENDDO
278    
279          DO k=1,Nr  C--     Start computation of dynamics
280           DO j=1-OLy,sNy+OLy          iMin = 1-OLx+2
281            DO i=1-OLx,sNx+OLx          iMax = sNx+OLx-1
282  C This is currently also used by IVDC and Diagnostics          jMin = 1-OLy+2
283             ConvectCount(i,j,k) = 0.          jMax = sNy+OLy-1
            KappaRT(i,j,k) = 0. _d 0  
            KappaRS(i,j,k) = 0. _d 0  
           ENDDO  
          ENDDO  
         ENDDO  
   
         iMin = 1-OLx+1  
         iMax = sNx+OLx  
         jMin = 1-OLy+1  
         jMax = sNy+OLy  
   
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
 CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
 CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
 CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
 #ifdef ALLOW_PASSIVE_TRACER  
 CADJ STORE tr1  (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
 #endif  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 C--     Start of diagnostic loop  
         DO k=Nr,1,-1  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 C? Patrick, is this formula correct now that we change the loop range?  
 C? Do we still need this?  
 cph kkey formula corrected.  
 cph Needed for rhok, rhokm1, in the case useGMREDI.  
          kkey = (ikey-1)*Nr + k  
 CADJ STORE rhokm1(:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte  
 CADJ STORE rhok  (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 C--       Integrate continuity vertically for vertical velocity  
           CALL INTEGRATE_FOR_W(  
      I                         bi, bj, k, uVel, vVel,  
      O                         wVel,  
      I                         myThid )  
   
 #ifdef    ALLOW_OBCS  
 #ifdef    ALLOW_NONHYDROSTATIC  
 C--       Apply OBC to W if in N-H mode  
           IF (useOBCS.AND.nonHydrostatic) THEN  
             CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )  
           ENDIF  
 #endif    /* ALLOW_NONHYDROSTATIC */  
 #endif    /* ALLOW_OBCS */  
   
 C--       Calculate gradients of potential density for isoneutral  
 C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)  
 c         IF ( k.GT.1 .AND. (useGMRedi.OR.ivdc_kappa.NE.0.) ) THEN  
           IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.) ) THEN  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  
 CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
             CALL FIND_RHO(  
      I        bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,  
      I        theta, salt,  
      O        rhoK,  
      I        myThid )  
             IF (k.GT.1) THEN  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  
 CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
              CALL FIND_RHO(  
      I        bi, bj, iMin, iMax, jMin, jMax, k-1, k, eosType,  
      I        theta, salt,  
      O        rhoKm1,  
      I        myThid )  
             ENDIF  
             CALL GRAD_SIGMA(  
      I             bi, bj, iMin, iMax, jMin, jMax, k,  
      I             rhoK, rhoKm1, rhoK,  
      O             sigmaX, sigmaY, sigmaR,  
      I             myThid )  
           ENDIF  
   
 C--       Implicit Vertical Diffusion for Convection  
 c ==> should use sigmaR !!!  
           IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN  
             CALL CALC_IVDC(  
      I        bi, bj, iMin, iMax, jMin, jMax, k,  
      I        rhoKm1, rhoK,  
      U        ConvectCount, KappaRT, KappaRS,  
      I        myTime, myIter, myThid)  
           ENDIF  
   
 C--     end of diagnostic k loop (Nr:1)  
         ENDDO  
284    
285  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
 cph avoids recomputation of integrate_for_w  
286  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
287  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
288    
289  #ifdef  ALLOW_OBCS  C--     Explicit part of the Surface Potentiel Gradient (add in TIMESTEP)
290  C--     Calculate future values on open boundaries  C       (note: this loop will be replaced by CALL CALC_GRAD_ETA)
291          IF (useOBCS) THEN          IF (implicSurfPress.NE.1.) THEN
292            CALL OBCS_CALC( bi, bj, myTime+deltaT,            CALL CALC_GRAD_PHI_SURF(
293       I            uVel, vVel, wVel, theta, salt,       I         bi,bj,iMin,iMax,jMin,jMax,
294       I            myThid )       I         etaN,
295          ENDIF       O         phiSurfX,phiSurfY,
296  #endif  /* ALLOW_OBCS */       I         myThid )                        
   
 C--     Determines forcing terms based on external fields  
 C       relaxation terms, etc.  
         CALL EXTERNAL_FORCING_SURF(  
      I             bi, bj, iMin, iMax, jMin, jMax,  
      I             myThid )  
 #ifdef ALLOW_AUTODIFF_TAMC  
 cph needed for KPP  
 CADJ STORE surfacetendencyU(:,:,bi,bj)  
 CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE surfacetendencyV(:,:,bi,bj)  
 CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE surfacetendencyS(:,:,bi,bj)  
 CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE surfacetendencyT(:,:,bi,bj)  
 CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 #ifdef  ALLOW_GMREDI  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE sigmaX(:,:,:) = comlev1, key=ikey, byte=isbyte  
 CADJ STORE sigmaY(:,:,:) = comlev1, key=ikey, byte=isbyte  
 CADJ STORE sigmaR(:,:,:) = comlev1, key=ikey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
 C--     Calculate iso-neutral slopes for the GM/Redi parameterisation  
         IF (useGMRedi) THEN  
           DO k=1,Nr  
             CALL GMREDI_CALC_TENSOR(  
      I             bi, bj, iMin, iMax, jMin, jMax, k,  
      I             sigmaX, sigmaY, sigmaR,  
      I             myThid )  
           ENDDO  
 #ifdef ALLOW_AUTODIFF_TAMC  
         ELSE  
           DO k=1, Nr  
             CALL GMREDI_CALC_TENSOR_DUMMY(  
      I             bi, bj, iMin, iMax, jMin, jMax, k,  
      I             sigmaX, sigmaY, sigmaR,  
      I             myThid )  
           ENDDO  
 #endif /* ALLOW_AUTODIFF_TAMC */  
         ENDIF  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE Kwx(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE Kwy(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE Kwz(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 #endif  /* ALLOW_GMREDI */  
   
 #ifdef  ALLOW_KPP  
 C--     Compute KPP mixing coefficients  
         IF (useKPP) THEN  
           CALL KPP_CALC(  
      I                  bi, bj, myTime, myThid )  
 #ifdef ALLOW_AUTODIFF_TAMC  
         ELSE  
           CALL KPP_CALC_DUMMY(  
      I                  bi, bj, myTime, myThid )  
 #endif /* ALLOW_AUTODIFF_TAMC */  
297          ENDIF          ENDIF
298    
299  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
300  CADJ STORE KPPghat   (:,:,:,bi,bj)  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
301  CADJ &   , KPPviscAz (:,:,:,bi,bj)  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
302  CADJ &   , KPPdiffKzT(:,:,:,bi,bj)  #ifdef ALLOW_KPP
303  CADJ &   , KPPdiffKzS(:,:,:,bi,bj)  CADJ STORE KPPviscAz (:,:,:,bi,bj)
 CADJ &   , KPPfrac   (:,:  ,bi,bj)  
304  CADJ &                 = comlev1_bibj, key=ikey, byte=isbyte  CADJ &                 = comlev1_bibj, key=ikey, byte=isbyte
305  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_KPP */
   
 #endif  /* ALLOW_KPP */  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE KappaRT(:,:,:)     = comlev1_bibj, key = ikey, byte = isbyte  
 CADJ STORE KappaRS(:,:,:)     = comlev1_bibj, key = ikey, byte = isbyte  
 CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
 CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
 CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
 CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
 #ifdef ALLOW_PASSIVE_TRACER  
 CADJ STORE tr1  (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
 #endif  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 #ifdef ALLOW_AIM  
 C       AIM - atmospheric intermediate model, physics package code.  
 C note(jmc) : phiHyd=0 at this point but is not really used in Molteni Physics  
         IF ( useAIM ) THEN  
          CALL TIMER_START('AIM_DO_ATMOS_PHYS      [DYNAMICS]', myThid)  
          CALL AIM_DO_ATMOS_PHYSICS( phiHyd, bi, bj, myTime, myThid )  
          CALL TIMER_STOP ('AIM_DO_ATMOS_PHYS      [DYNAMICS]', myThid)  
         ENDIF  
 #endif /* ALLOW_AIM */  
   
   
 C--     Start of thermodynamics loop  
         DO k=Nr,1,-1  
 #ifdef ALLOW_AUTODIFF_TAMC  
 C? Patrick Is this formula correct?  
 cph Yes, but I rewrote it.  
 cph Also, the KappaR? need the index and subscript k!  
          kkey = (ikey-1)*Nr + k  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 C--       km1    Points to level above k (=k-1)  
 C--       kup    Cycles through 1,2 to point to layer above  
 C--       kDown  Cycles through 2,1 to point to current layer  
   
           km1  = MAX(1,k-1)  
           kup  = 1+MOD(k+1,2)  
           kDown= 1+MOD(k,2)  
   
           iMin = 1-OLx  
           iMax = sNx+OLx  
           jMin = 1-OLy  
           jMax = sNy+OLy  
   
 C--       Get temporary terms used by tendency routines  
           CALL CALC_COMMON_FACTORS (  
      I         bi,bj,iMin,iMax,jMin,jMax,k,  
      O         xA,yA,uTrans,vTrans,rTrans,maskUp,  
      I         myThid)  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE KappaRT(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte  
 CADJ STORE KappaRS(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte  
306  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
307    
308  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL
309  C--      Calculate the total vertical diffusivity  C--      Calculate the total vertical diffusivity
310           CALL CALC_DIFFUSIVITY(          DO k=1,Nr
311             CALL CALC_VISCOSITY(
312       I        bi,bj,iMin,iMax,jMin,jMax,k,       I        bi,bj,iMin,iMax,jMin,jMax,k,
313       I        maskUp,       O        KappaRU,KappaRV,
      O        KappaRT,KappaRS,KappaRU,KappaRV,  
314       I        myThid)       I        myThid)
315           ENDDO
316  #endif  #endif
317    
           iMin = 1-OLx+2  
           iMax = sNx+OLx-1  
           jMin = 1-OLy+2  
           jMax = sNy+OLy-1  
   
 C--      Calculate active tracer tendencies (gT,gS,...)  
 C        and step forward storing result in gTnm1, gSnm1, etc.  
          IF ( tempStepping ) THEN  
            CALL CALC_GT(  
      I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,  
      I         xA,yA,uTrans,vTrans,rTrans,maskUp,  
      I         KappaRT,  
      U         fVerT,  
      I         myTime, myThid)  
            tauAB = 0.5d0 + abEps  
            CALL TIMESTEP_TRACER(  
      I         bi,bj,iMin,iMax,jMin,jMax,k,tauAB,  
      I         theta, gT,  
      U         gTnm1,  
      I         myIter, myThid)  
          ENDIF  
          IF ( saltStepping ) THEN  
            CALL CALC_GS(  
      I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,  
      I         xA,yA,uTrans,vTrans,rTrans,maskUp,  
      I         KappaRS,  
      U         fVerS,  
      I         myTime, myThid)  
            tauAB = 0.5d0 + abEps  
            CALL TIMESTEP_TRACER(  
      I         bi,bj,iMin,iMax,jMin,jMax,k,tauAB,  
      I         salt, gS,  
      U         gSnm1,  
      I         myIter, myThid)  
          ENDIF  
 #ifdef ALLOW_PASSIVE_TRACER  
          IF ( tr1Stepping ) THEN  
            CALL CALC_GTR1(  
      I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,  
      I         xA,yA,uTrans,vTrans,rTrans,maskUp,  
      I         KappaRT,  
      U         fVerTr1,  
      I         myTime, myThid)  
            tauAB = 0.5d0 + abEps  
            CALL TIMESTEP_TRACER(  
      I         bi,bj,iMin,iMax,jMin,jMax,k,tauAB,  
      I         Tr1, gTr1,  
      U         gTr1NM1,  
      I         myIter, myThid)  
          ENDIF  
 #endif  
   
 #ifdef   ALLOW_OBCS  
 C--      Apply open boundary conditions  
          IF (useOBCS) THEN  
            CALL OBCS_APPLY_TS( bi, bj, k, gTnm1, gSnm1, myThid )  
          END IF  
 #endif   /* ALLOW_OBCS */  
   
 C--      Freeze water  
          IF (allowFreezing) THEN  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE gTNm1(:,:,k,bi,bj) = comlev1_bibj_k  
 CADJ &   , key = kkey, byte = isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
             CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )  
          END IF  
   
 C--     end of thermodynamic k loop (Nr:1)  
         ENDDO  
   
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 C? Patrick? What about this one?  
 cph Keys iikey and idkey don't seem to be needed  
 cph since storing occurs on different tape for each  
 cph impldiff call anyways.  
 cph Thus, common block comlev1_impl isn't needed either.  
 cph Storing below needed in the case useGMREDI.  
         iikey = (ikey-1)*maximpl  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 C--     Implicit diffusion  
         IF (implicitDiffusion) THEN  
   
          IF (tempStepping) THEN  
 #ifdef ALLOW_AUTODIFF_TAMC  
             idkey = iikey + 1  
 CADJ STORE gTNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
             CALL IMPLDIFF(  
      I         bi, bj, iMin, iMax, jMin, jMax,  
      I         deltaTtracer, KappaRT, recip_HFacC,  
      U         gTNm1,  
      I         myThid )  
          ENDIF  
   
          IF (saltStepping) THEN  
 #ifdef ALLOW_AUTODIFF_TAMC  
          idkey = iikey + 2  
 CADJ STORE gSNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
             CALL IMPLDIFF(  
      I         bi, bj, iMin, iMax, jMin, jMax,  
      I         deltaTtracer, KappaRS, recip_HFacC,  
      U         gSNm1,  
      I         myThid )  
          ENDIF  
   
 #ifdef ALLOW_PASSIVE_TRACER  
          IF (tr1Stepping) THEN  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE gTr1Nm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
           CALL IMPLDIFF(  
      I      bi, bj, iMin, iMax, jMin, jMax,  
      I      deltaTtracer, KappaRT, recip_HFacC,  
      U      gTr1Nm1,  
      I      myThid )  
          ENDIF  
 #endif  
   
 #ifdef   ALLOW_OBCS  
 C--      Apply open boundary conditions  
          IF (useOBCS) THEN  
            DO K=1,Nr  
              CALL OBCS_APPLY_TS( bi, bj, k, gTnm1, gSnm1, myThid )  
            ENDDO  
          END IF  
 #endif   /* ALLOW_OBCS */  
   
 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 Potentiel Gradient (add in TIMESTEP)  
 C       (note: this loop will be replaced by CALL CALC_GRAD_ETA)  
         IF (implicSurfPress.NE.1.) THEN  
           CALL CALC_GRAD_PHI_SURF(  
      I         bi,bj,iMin,iMax,jMin,jMax,  
      I         etaN,  
      O         phiSurfX,phiSurfY,  
      I         myThid )                          
         ENDIF  
   
318  C--     Start of dynamics loop  C--     Start of dynamics loop
319          DO k=1,Nr          DO k=1,Nr
320    
# Line 663  C--       kup    Cycles through 1,2 to p Line 323  C--       kup    Cycles through 1,2 to p
323  C--       kDown  Cycles through 2,1 to point to current layer  C--       kDown  Cycles through 2,1 to point to current layer
324    
325            km1  = MAX(1,k-1)            km1  = MAX(1,k-1)
326              kp1  = MIN(k+1,Nr)
327            kup  = 1+MOD(k+1,2)            kup  = 1+MOD(k+1,2)
328            kDown= 1+MOD(k,2)            kDown= 1+MOD(k,2)
329    
330    #ifdef ALLOW_AUTODIFF_TAMC
331             kkey = (ikey-1)*Nr + k
332    #endif /* ALLOW_AUTODIFF_TAMC */
333    
334  C--      Integrate hydrostatic balance for phiHyd with BC of  C--      Integrate hydrostatic balance for phiHyd with BC of
335  C        phiHyd(z=0)=0  C        phiHyd(z=0)=0
336  C        distinguishe between Stagger and Non Stagger time stepping  C        distinguishe between Stagger and Non Stagger time stepping
337           IF (staggerTimeStep) THEN           IF (staggerTimeStep) THEN
338             CALL CALC_PHI_HYD(             CALL CALC_PHI_HYD(
339       I        bi,bj,iMin,iMax,jMin,jMax,k,       I        bi,bj,iMin,iMax,jMin,jMax,k,
340       I        gTnm1, gSnm1,       I        gT, gS,
341       U        phiHyd,       U        phiHyd,
342       I        myThid )       I        myThid )
343           ELSE           ELSE
# Line 683  C        distinguishe between Stagger an Line 348  C        distinguishe between Stagger an
348       I        myThid )       I        myThid )
349           ENDIF           ENDIF
350    
351    C        calculate pressure from phiHyd and store it on common block
352    C        variable pressure
353             CALL STORE_PRESSURE( bi, bj, k, phiHyd, myThid )
354    
355    
356  C--      Calculate accelerations in the momentum equations (gU, gV, ...)  C--      Calculate accelerations in the momentum equations (gU, gV, ...)
357  C        and step forward storing the result in gUnm1, gVnm1, etc...  C        and step forward storing the result in gUnm1, gVnm1, etc...
358           IF ( momStepping ) THEN           IF ( momStepping ) THEN
359             CALL CALC_MOM_RHS(  #ifndef DISABLE_MOM_FLUXFORM
360               IF (.NOT. vectorInvariantMomentum) CALL MOM_FLUXFORM(
361         I         bi,bj,iMin,iMax,jMin,jMax,k,kup,kDown,
362         I         phiHyd,KappaRU,KappaRV,
363         U         fVerU, fVerV,
364         I         myTime, myIter, myThid)
365    #endif
366    #ifndef DISABLE_MOM_VECINV
367               IF (vectorInvariantMomentum) CALL MOM_VECINV(
368       I         bi,bj,iMin,iMax,jMin,jMax,k,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax,k,kup,kDown,
369       I         phiHyd,KappaRU,KappaRV,       I         phiHyd,KappaRU,KappaRV,
370       U         fVerU, fVerV,       U         fVerU, fVerV,
371       I         myTime, myThid)       I         myTime, myIter, myThid)
372    #endif
373             CALL TIMESTEP(             CALL TIMESTEP(
374       I         bi,bj,iMin,iMax,jMin,jMax,k,       I         bi,bj,iMin,iMax,jMin,jMax,k,
375       I         phiHyd, phiSurfX, phiSurfY,       I         phiHyd, phiSurfX, phiSurfY,
# Line 720  C--      Apply open boundary conditions Line 399  C--      Apply open boundary conditions
399  C--     end of dynamics k loop (1:Nr)  C--     end of dynamics k loop (1:Nr)
400          ENDDO          ENDDO
401    
   
   
402  C--     Implicit viscosity  C--     Implicit viscosity
403          IF (implicitViscosity.AND.momStepping) THEN          IF (implicitViscosity.AND.momStepping) THEN
404  #ifdef    ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
           idkey = iikey + 3  
405  CADJ STORE gUNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE gUNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
406  #endif    /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
407            CALL IMPLDIFF(            CALL IMPLDIFF(
# Line 734  CADJ STORE gUNm1(:,:,:,bi,bj) = comlev1_ Line 410  CADJ STORE gUNm1(:,:,:,bi,bj) = comlev1_
410       U         gUNm1,       U         gUNm1,
411       I         myThid )       I         myThid )
412  #ifdef    ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
           idkey = iikey + 4  
413  CADJ STORE gVNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE gVNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
414  #endif    /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
415            CALL IMPLDIFF(            CALL IMPLDIFF(
# Line 754  C--      Apply open boundary conditions Line 429  C--      Apply open boundary conditions
429    
430  #ifdef    INCLUDE_CD_CODE  #ifdef    INCLUDE_CD_CODE
431  #ifdef    ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
           idkey = iikey + 5  
432  CADJ STORE vVelD(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE vVelD(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
433  #endif    /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
434            CALL IMPLDIFF(            CALL IMPLDIFF(
# Line 763  CADJ STORE vVelD(:,:,:,bi,bj) = comlev1_ Line 437  CADJ STORE vVelD(:,:,:,bi,bj) = comlev1_
437       U         vVelD,       U         vVelD,
438       I         myThid )       I         myThid )
439  #ifdef    ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
           idkey = iikey + 6  
440  CADJ STORE uVelD(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE uVelD(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
441  #endif    /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
442            CALL IMPLDIFF(            CALL IMPLDIFF(
# Line 776  C--     End If implicitViscosity.AND.mom Line 449  C--     End If implicitViscosity.AND.mom
449          ENDIF          ENDIF
450    
451  Cjmc : add for phiHyd output <- but not working if multi tile per CPU  Cjmc : add for phiHyd output <- but not working if multi tile per CPU
452  c       IF ( DIFFERENT_MULTIPLE(dumpFreq,myTime+deltaTClock,myTime)  c        IF ( DIFFERENT_MULTIPLE(dumpFreq,myTime+deltaTClock,myTime)
453  c    &  .AND. buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN  c     &       .AND. buoyancyRelation .ne. 'OCEANIC' ) THEN
454  c         WRITE(suff,'(I10.10)') myIter+1  c           WRITE(suff,'(I10.10)') myIter+1
455  c         CALL WRITE_FLD_XYZ_RL('PH.',suff,phiHyd,myIter+1,myThid)  c           CALL WRITE_FLD_XYZ_RL('PH.',suff,phiHyd,myIter+1,myThid)
456  c       ENDIF  c        ENDIF
457  Cjmc(end)  Cjmc(end)
458    
459  #ifdef ALLOW_TIMEAVE  #ifdef ALLOW_TIMEAVE
460          IF (taveFreq.GT.0.) THEN          IF (taveFreq.GT.0.) THEN
461            CALL TIMEAVE_CUMUL_1T(phiHydtave, phiHyd, Nr,            CALL TIMEAVE_CUMUL_1T(phiHydtave, phiHyd, Nr,
462       I                              deltaTclock, bi, bj, myThid)       I                              deltaTclock, bi, bj, myThid)
           IF (ivdc_kappa.NE.0.) THEN  
             CALL TIMEAVE_CUMULATE(ConvectCountTave, ConvectCount, Nr,  
      I                              deltaTclock, bi, bj, myThid)  
           ENDIF  
463          ENDIF          ENDIF
464  #endif /* ALLOW_TIMEAVE */  #endif /* ALLOW_TIMEAVE */
465    
466         ENDDO         ENDDO
467        ENDDO        ENDDO
468    
469  #ifndef EXCLUDE_DEBUGMODE  Cml(
470    C     In order to compare the variance of phiHydLow of a p/z-coordinate
471    C     run with etaH of a z/p-coordinate run the drift of phiHydLow
472    C     has to be removed by something like the following subroutine:
473    C      CALL REMOVE_MEAN_RL( 1, phiHydLow, maskH, maskH, rA, drF,
474    C     &                'phiHydLow', myThid )
475    Cml)
476    
477    #ifndef DISABLE_DEBUGMODE
478        If (debugMode) THEN        If (debugMode) THEN
479         CALL DEBUG_STATS_RL(1,EtaN,'EtaN (DYNAMICS)',myThid)         CALL DEBUG_STATS_RL(1,EtaN,'EtaN (DYNAMICS)',myThid)
480         CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (DYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (DYNAMICS)',myThid)
Legend:
Removed from v.1.74  
changed lines
  Added in v.1.90
  ViewVC Help
Powered by ViewVC 1.1.22