/[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.68 by adcroft, Tue May 29 14:01:37 2001 UTC revision 1.103 by edhill, Thu Oct 30 12:00:41 2003 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    
4    #include "PACKAGES_CONFIG.h"
5  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
6    
7    CBOP
8    C     !ROUTINE: DYNAMICS
9    C     !INTERFACE:
10        SUBROUTINE DYNAMICS(myTime, myIter, myThid)        SUBROUTINE DYNAMICS(myTime, myIter, myThid)
11  C     /==========================================================\  C     !DESCRIPTION: \bv
12  C     | SUBROUTINE DYNAMICS                                      |  C     *==========================================================*
13  C     | o Controlling routine for the explicit part of the model |  C     | SUBROUTINE DYNAMICS                                      
14  C     |   dynamics.                                              |  C     | o Controlling routine for the explicit part of the model  
15  C     |==========================================================|  C     |   dynamics.                                              
16  C     | This routine evaluates the "dynamics" terms for each     |  C     *==========================================================*
17  C     | block of ocean in turn. Because the blocks of ocean have |  C     | This routine evaluates the "dynamics" terms for each      
18  C     | overlap regions they are independent of one another.     |  C     | block of ocean in turn. Because the blocks of ocean have  
19  C     | If terms involving lateral integrals are needed in this  |  C     | overlap regions they are independent of one another.      
20  C     | routine care will be needed. Similarly finite-difference |  C     | If terms involving lateral integrals are needed in this  
21  C     | operations with stencils wider than the overlap region   |  C     | routine care will be needed. Similarly finite-difference  
22  C     | require special consideration.                           |  C     | operations with stencils wider than the overlap region    
23  C     | Notes                                                    |  C     | require special consideration.                            
24  C     | =====                                                    |  C     | The algorithm...
25  C     | C*P* comments indicating place holders for which code is |  C     |
26  C     |      presently being developed.                          |  C     | "Correction Step"
27  C     \==========================================================/  C     | =================
28    C     | Here we update the horizontal velocities with the surface
29    C     | pressure such that the resulting flow is either consistent
30    C     | with the free-surface evolution or the rigid-lid:
31    C     |   U[n] = U* + dt x d/dx P
32    C     |   V[n] = V* + dt x d/dy P
33    C     |
34    C     | "Calculation of Gs"
35    C     | ===================
36    C     | This is where all the accelerations and tendencies (ie.
37    C     | physics, parameterizations etc...) are calculated
38    C     |   rho = rho ( theta[n], salt[n] )
39    C     |   b   = b(rho, theta)
40    C     |   K31 = K31 ( rho )
41    C     |   Gu[n] = Gu( u[n], v[n], wVel, b, ... )
42    C     |   Gv[n] = Gv( u[n], v[n], wVel, b, ... )
43    C     |   Gt[n] = Gt( theta[n], u[n], v[n], wVel, K31, ... )
44    C     |   Gs[n] = Gs( salt[n], u[n], v[n], wVel, K31, ... )
45    C     |
46    C     | "Time-stepping" or "Prediction"
47    C     | ================================
48    C     | The models variables are stepped forward with the appropriate
49    C     | time-stepping scheme (currently we use Adams-Bashforth II)
50    C     | - For momentum, the result is always *only* a "prediction"
51    C     | in that the flow may be divergent and will be "corrected"
52    C     | later with a surface pressure gradient.
53    C     | - Normally for tracers the result is the new field at time
54    C     | level [n+1} *BUT* in the case of implicit diffusion the result
55    C     | is also *only* a prediction.
56    C     | - We denote "predictors" with an asterisk (*).
57    C     |   U* = U[n] + dt x ( 3/2 Gu[n] - 1/2 Gu[n-1] )
58    C     |   V* = V[n] + dt x ( 3/2 Gv[n] - 1/2 Gv[n-1] )
59    C     |   theta[n+1] = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
60    C     |   salt[n+1] = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
61    C     | With implicit diffusion:
62    C     |   theta* = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
63    C     |   salt* = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
64    C     |   (1 + dt * K * d_zz) theta[n] = theta*
65    C     |   (1 + dt * K * d_zz) salt[n] = salt*
66    C     |
67    C     *==========================================================*
68    C     \ev
69    C     !USES:
70        IMPLICIT NONE        IMPLICIT NONE
   
71  C     == Global variables ===  C     == Global variables ===
72  #include "SIZE.h"  #include "SIZE.h"
73  #include "EEPARAMS.h"  #include "EEPARAMS.h"
74  #include "PARAMS.h"  #include "PARAMS.h"
75  #include "DYNVARS.h"  #include "DYNVARS.h"
76    #ifdef ALLOW_CD_CODE
77    #include "CD_CODE_VARS.h"
78    #endif
79  #include "GRID.h"  #include "GRID.h"
80    #ifdef ALLOW_PASSIVE_TRACER
81    #include "TR1.h"
82    #endif
83  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
84  # include "tamc.h"  # include "tamc.h"
85  # include "tamc_keys.h"  # include "tamc_keys.h"
86  # include "FFIELDS.h"  # include "FFIELDS.h"
87    # include "EOS.h"
88  # ifdef ALLOW_KPP  # ifdef ALLOW_KPP
89  #  include "KPP.h"  #  include "KPP.h"
90  # endif  # endif
 # ifdef ALLOW_GMREDI  
 #  include "GMREDI.h"  
 # endif  
91  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
92    
93  #ifdef ALLOW_TIMEAVE  C     !CALLING SEQUENCE:
94  #include "TIMEAVE_STATV.h"  C     DYNAMICS()
95  #endif  C      |
96    C      |-- CALC_GRAD_PHI_SURF
97    C      |
98    C      |-- CALC_VISCOSITY
99    C      |
100    C      |-- CALC_PHI_HYD  
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 54  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
131  C                                      variable.  C                                      variable.
132  C     rhoK, rhoKM1   - Density at current level, and level above  C     phiHydC    :: hydrostatic potential anomaly at cell center
133  C     phiHyd         - Hydrostatic part of the potential phiHydi.  C                   In z coords phiHyd is the hydrostatic potential
134  C                      In z coords phiHydiHyd is the hydrostatic  C                      (=pressure/rho0) anomaly
135  C                      Potential (=pressure/rho0) anomaly  C                   In p coords phiHyd is the geopotential height anomaly.
136  C                      In p coords phiHydiHyd is the geopotential  C     phiHydF    :: hydrostatic potential anomaly at middle between 2 centers
137  C                      surface height anomaly.  C     dPhiHydX,Y :: Gradient (X & Y directions) of hydrostatic potential anom.
138  C     phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean)  C     phiSurfX,  ::  gradient of Surface potential (Pressure/rho, ocean)
139  C     phiSurfY             or geopotentiel (atmos) in X and Y direction  C     phiSurfY             or geopotential (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)  
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 phiHydF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
149        _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiHydC (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
150        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL dPhiHydX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
151          _RL dPhiHydY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
152        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
153        _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)  
154        _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
155        _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)  
156    
157        INTEGER iMin, iMax        INTEGER iMin, iMax
158        INTEGER jMin, jMax        INTEGER jMin, jMax
159        INTEGER bi, bj        INTEGER bi, bj
160        INTEGER i, j        INTEGER i, j
161        INTEGER k, km1, kup, kDown        INTEGER k, km1, kp1, kup, kDown
162    
163  Cjmc : add for phiHyd output <- but not working if multi tile per CPU        LOGICAL  DIFFERENT_MULTIPLE
164  c     CHARACTER*(MAX_LEN_MBUF) suff        EXTERNAL DIFFERENT_MULTIPLE
 c     LOGICAL  DIFFERENT_MULTIPLE  
 c     EXTERNAL DIFFERENT_MULTIPLE  
 Cjmc(end)  
165    
166  C---    The algorithm...  C---    The algorithm...
167  C  C
# Line 165  C         salt* = salt[n] + dt x ( 3/2 G Line 206  C         salt* = salt[n] + dt x ( 3/2 G
206  C         (1 + dt * K * d_zz) theta[n] = theta*  C         (1 + dt * K * d_zz) theta[n] = theta*
207  C         (1 + dt * K * d_zz) salt[n] = salt*  C         (1 + dt * K * d_zz) salt[n] = salt*
208  C---  C---
209    CEOP
210    
211  #ifdef ALLOW_AUTODIFF_TAMC  C-- Call to routine for calculation of
212  C--   dummy statement to end declaration part  C   Eliassen-Palm-flux-forced U-tendency,
213        ikey = 1  C   if desired:
214  #endif /* ALLOW_AUTODIFF_TAMC */  #ifdef INCLUDE_EP_FORCING_CODE
215          CALL CALC_EP_FORCING(myThid)
216  C--   Set up work arrays with valid (i.e. not NaN) values  #endif
 C     These inital values do not alter the numerical results. They  
 C     just ensure that all memory references are to valid floating  
 C     point numbers. This prevents spurious hardware signals due to  
 C     uninitialised but inert locations.  
       DO j=1-OLy,sNy+OLy  
        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  
         rhoKM1 (i,j) = 0. _d 0  
         rhok   (i,j) = 0. _d 0  
         phiSurfX(i,j) = 0. _d 0  
         phiSurfY(i,j) = 0. _d 0  
        ENDDO  
       ENDDO  
   
217    
218  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
219  C--   HPF directive to help TAMC  C--   HPF directive to help TAMC
# Line 207  CHPF$ INDEPENDENT Line 224  CHPF$ INDEPENDENT
224    
225  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
226  C--    HPF directive to help TAMC  C--    HPF directive to help TAMC
227  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS,fVerU,fVerV  CHPF$  INDEPENDENT, NEW (fVerU,fVerV
228  CHPF$&                  ,phiHyd,utrans,vtrans,xA,yA  CHPF$&                  ,phiHydF
229  CHPF$&                  ,KappaRT,KappaRS,KappaRU,KappaRV  CHPF$&                  ,KappaRU,KappaRV
230  CHPF$&                  )  CHPF$&                  )
231  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
232    
# Line 218  CHPF$&                  ) Line 235  CHPF$&                  )
235  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
236            act1 = bi - myBxLo(myThid)            act1 = bi - myBxLo(myThid)
237            max1 = myBxHi(myThid) - myBxLo(myThid) + 1            max1 = myBxHi(myThid) - myBxLo(myThid) + 1
   
238            act2 = bj - myByLo(myThid)            act2 = bj - myByLo(myThid)
239            max2 = myByHi(myThid) - myByLo(myThid) + 1            max2 = myByHi(myThid) - myByLo(myThid) + 1
   
240            act3 = myThid - 1            act3 = myThid - 1
241            max3 = nTx*nTy            max3 = nTx*nTy
   
242            act4 = ikey_dynamics - 1            act4 = ikey_dynamics - 1
243              idynkey = (act1 + 1) + act2*max1
           ikey = (act1 + 1) + act2*max1  
244       &                      + act3*max1*max2       &                      + act3*max1*max2
245       &                      + act4*max1*max2*max3       &                      + act4*max1*max2*max3
246  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
247    
248  C--     Set up work arrays that need valid initial values  C--   Set up work arrays with valid (i.e. not NaN) values
249          DO j=1-OLy,sNy+OLy  C     These inital values do not alter the numerical results. They
250           DO i=1-OLx,sNx+OLx  C     just ensure that all memory references are to valid floating
251            rTrans(i,j)   = 0. _d 0  C     point numbers. This prevents spurious hardware signals due to
252            fVerT (i,j,1) = 0. _d 0  C     uninitialised but inert locations.
           fVerT (i,j,2) = 0. _d 0  
           fVerS (i,j,1) = 0. _d 0  
           fVerS (i,j,2) = 0. _d 0  
           fVerU (i,j,1) = 0. _d 0  
           fVerU (i,j,2) = 0. _d 0  
           fVerV (i,j,1) = 0. _d 0  
           fVerV (i,j,2) = 0. _d 0  
          ENDDO  
         ENDDO  
253    
254          DO k=1,Nr          DO k=1,Nr
255           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
256            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
257  C This is currently also used by IVDC and Diagnostics             KappaRU(i,j,k) = 0. _d 0
258             ConvectCount(i,j,k) = 0.             KappaRV(i,j,k) = 0. _d 0
259             KappaRT(i,j,k) = 0. _d 0  #ifdef ALLOW_AUTODIFF_TAMC
260             KappaRS(i,j,k) = 0. _d 0  cph(
261    c--   need some re-initialisation here to break dependencies
262    c--   totphihyd is assumed zero from ini_pressure, i.e.
263    c--   avoiding iterate pressure p = integral of (g*rho(p)*dz)
264    cph)
265               totPhiHyd(i,j,k,bi,bj) = 0. _d 0
266               gu(i,j,k,bi,bj) = 0. _d 0
267               gv(i,j,k,bi,bj) = 0. _d 0
268    #endif
269            ENDDO            ENDDO
270           ENDDO           ENDDO
271          ENDDO          ENDDO
272            DO j=1-OLy,sNy+OLy
273          iMin = 1-OLx+1           DO i=1-OLx,sNx+OLx
274          iMax = sNx+OLx            fVerU  (i,j,1) = 0. _d 0
275          jMin = 1-OLy+1            fVerU  (i,j,2) = 0. _d 0
276          jMax = sNy+OLy            fVerV  (i,j,1) = 0. _d 0
277              fVerV  (i,j,2) = 0. _d 0
278              phiHydF (i,j)  = 0. _d 0
279  #ifdef ALLOW_AUTODIFF_TAMC            phiHydC (i,j)  = 0. _d 0
280  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte            dPhiHydX(i,j)  = 0. _d 0
281  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte            dPhiHydY(i,j)  = 0. _d 0
282  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte            phiSurfX(i,j)  = 0. _d 0
283  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte            phiSurfY(i,j)  = 0. _d 0
284  #endif /* ALLOW_AUTODIFF_TAMC */           ENDDO
   
 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)  
285          ENDDO          ENDDO
286    
287  #ifdef ALLOW_AUTODIFF_TAMC  C--     Start computation of dynamics
288  cph avoids recomputation of integrate_for_w          iMin = 0
289  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte          iMax = sNx+1
290  #endif /* ALLOW_AUTODIFF_TAMC */          jMin = 0
291            jMax = sNy+1
 #ifdef  ALLOW_OBCS  
 C--     Calculate future values on open boundaries  
         IF (useOBCS) THEN  
           CALL OBCS_CALC( bi, bj, myTime+deltaT,  
      I            uVel, vVel, wVel, theta, salt,  
      I            myThid )  
         ENDIF  
 #endif  /* ALLOW_OBCS */  
   
 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  
292    
293  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
294  CADJ STORE Kwx(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE wvel (:,:,:,bi,bj) =
295  CADJ STORE Kwy(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte  CADJ &     comlev1_bibj, key = idynkey, byte = isbyte
 CADJ STORE Kwz(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte  
296  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
297    
298  #endif  /* ALLOW_GMREDI */  C--     Explicit part of the Surface Potentiel Gradient (add in TIMESTEP)
299    C       (note: this loop will be replaced by CALL CALC_GRAD_ETA)
300  #ifdef  ALLOW_KPP          IF (implicSurfPress.NE.1.) THEN
301  C--     Compute KPP mixing coefficients            CALL CALC_GRAD_PHI_SURF(
302          IF (useKPP) THEN       I         bi,bj,iMin,iMax,jMin,jMax,
303            CALL KPP_CALC(       I         etaN,
304       I                  bi, bj, myTime, myThid )       O         phiSurfX,phiSurfY,
305  #ifdef ALLOW_AUTODIFF_TAMC       I         myThid )                        
         ELSE  
           CALL KPP_CALC_DUMMY(  
      I                  bi, bj, myTime, myThid )  
 #endif /* ALLOW_AUTODIFF_TAMC */  
306          ENDIF          ENDIF
307    
308  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
309  CADJ STORE KPPghat   (:,:,:,bi,bj)  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=idynkey, byte=isbyte
310  CADJ &   , KPPviscAz (:,:,:,bi,bj)  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=idynkey, byte=isbyte
311  CADJ &   , KPPdiffKzT(:,:,:,bi,bj)  #ifdef ALLOW_KPP
312  CADJ &   , KPPdiffKzS(:,:,:,bi,bj)  CADJ STORE KPPviscAz (:,:,:,bi,bj)
313  CADJ &   , KPPfrac   (:,:  ,bi,bj)  CADJ &                 = comlev1_bibj, key=idynkey, byte=isbyte
314  CADJ &                 = comlev1_bibj, key=ikey, byte=isbyte  #endif /* ALLOW_KPP */
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 #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  
 #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, 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+2  
           iMax = sNx+OLx-1  
           jMin = 1-OLy+2  
           jMax = sNy+OLy-1  
   
 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  
315  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
316    
317  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL
318  C--      Calculate the total vertical diffusivity  C--      Calculate the total vertical diffusivity
319           CALL CALC_DIFFUSIVITY(          DO k=1,Nr
320             CALL CALC_VISCOSITY(
321       I        bi,bj,iMin,iMax,jMin,jMax,k,       I        bi,bj,iMin,iMax,jMin,jMax,k,
322       I        maskUp,       O        KappaRU,KappaRV,
      O        KappaRT,KappaRS,KappaRU,KappaRV,  
323       I        myThid)       I        myThid)
324           ENDDO
325  #endif  #endif
326    
 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_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  
   
   
327  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
328  C? Patrick? What about this one?  CADJ STORE KappaRU(:,:,:)
329  cph Keys iikey and idkey don't seem to be needed  CADJ &                 = comlev1_bibj, key=idynkey, byte=isbyte
330  cph since storing occurs on different tape for each  CADJ STORE KappaRV(:,:,:)
331  cph impldiff call anyways.  CADJ &                 = comlev1_bibj, key=idynkey, byte=isbyte
 cph Thus, common block comlev1_impl isn't needed either.  
 cph Storing below needed in the case useGMREDI.  
         iikey = (ikey-1)*maximpl  
332  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
333    
 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_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  
   
334  C--     Start of dynamics loop  C--     Start of dynamics loop
335          DO k=1,Nr          DO k=1,Nr
336    
# Line 617  C--       kup    Cycles through 1,2 to p Line 339  C--       kup    Cycles through 1,2 to p
339  C--       kDown  Cycles through 2,1 to point to current layer  C--       kDown  Cycles through 2,1 to point to current layer
340    
341            km1  = MAX(1,k-1)            km1  = MAX(1,k-1)
342              kp1  = MIN(k+1,Nr)
343            kup  = 1+MOD(k+1,2)            kup  = 1+MOD(k+1,2)
344            kDown= 1+MOD(k,2)            kDown= 1+MOD(k,2)
345    
346    #ifdef ALLOW_AUTODIFF_TAMC
347             kkey = (idynkey-1)*Nr + k
348    c
349    CADJ STORE totphihyd (:,:,k,bi,bj)
350    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
351    CADJ STORE gt (:,:,k,bi,bj)
352    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
353    CADJ STORE gs (:,:,k,bi,bj)
354    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
355    CADJ STORE theta (:,:,k,bi,bj)
356    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
357    CADJ STORE salt  (:,:,k,bi,bj)
358    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
359    #endif /* ALLOW_AUTODIFF_TAMC */
360    
361  C--      Integrate hydrostatic balance for phiHyd with BC of  C--      Integrate hydrostatic balance for phiHyd with BC of
362  C        phiHyd(z=0)=0  C        phiHyd(z=0)=0
363  C        distinguishe between Stagger and Non Stagger time stepping  C        distinguishe between Stagger and Non Stagger time stepping
364           IF (staggerTimeStep) THEN           IF (staggerTimeStep) THEN
365             CALL CALC_PHI_HYD(             CALL CALC_PHI_HYD(
366       I        bi,bj,iMin,iMax,jMin,jMax,k,       I        bi,bj,iMin,iMax,jMin,jMax,k,
367       I        gTnm1, gSnm1,       I        gT, gS,
368       U        phiHyd,       U        phiHydF,
369       I        myThid )       O        phiHydC, dPhiHydX, dPhiHydY,
370         I        myTime, myIter, myThid )
371           ELSE           ELSE
372             CALL CALC_PHI_HYD(             CALL CALC_PHI_HYD(
373       I        bi,bj,iMin,iMax,jMin,jMax,k,       I        bi,bj,iMin,iMax,jMin,jMax,k,
374       I        theta, salt,       I        theta, salt,
375       U        phiHyd,       U        phiHydF,
376       I        myThid )       O        phiHydC, dPhiHydX, dPhiHydY,
377         I        myTime, myIter, myThid )
378           ENDIF           ENDIF
379    
380  C--      Calculate accelerations in the momentum equations (gU, gV, ...)  C--      Calculate accelerations in the momentum equations (gU, gV, ...)
381  C        and step forward storing the result in gUnm1, gVnm1, etc...  C        and step forward storing the result in gU, gV, etc...
382           IF ( momStepping ) THEN           IF ( momStepping ) THEN
383             CALL CALC_MOM_RHS(  #ifndef DISABLE_MOM_FLUXFORM
384               IF (.NOT. vectorInvariantMomentum) CALL MOM_FLUXFORM(
385         I         bi,bj,iMin,iMax,jMin,jMax,k,kup,kDown,
386         I         dPhiHydX,dPhiHydY,KappaRU,KappaRV,
387         U         fVerU, fVerV,
388         I         myTime, myIter, myThid)
389    #endif
390    #ifndef DISABLE_MOM_VECINV
391               IF (vectorInvariantMomentum) CALL MOM_VECINV(
392       I         bi,bj,iMin,iMax,jMin,jMax,k,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax,k,kup,kDown,
393       I         phiHyd,KappaRU,KappaRV,       I         dPhiHydX,dPhiHydY,KappaRU,KappaRV,
394       U         fVerU, fVerV,       U         fVerU, fVerV,
395       I         myTime, myThid)       I         myTime, myIter, myThid)
396    #endif
397             CALL TIMESTEP(             CALL TIMESTEP(
398       I         bi,bj,iMin,iMax,jMin,jMax,k,       I         bi,bj,iMin,iMax,jMin,jMax,k,
399       I         phiHyd, phiSurfX, phiSurfY,       I         dPhiHydX,dPhiHydY, phiSurfX, phiSurfY,
400       I         myIter, myThid)       I         myTime, myIter, myThid)
401    
402  #ifdef   ALLOW_OBCS  #ifdef   ALLOW_OBCS
403  C--      Apply open boundary conditions  C--      Apply open boundary conditions
404           IF (useOBCS) THEN             IF (useOBCS) THEN
405             CALL OBCS_APPLY_UV( bi, bj, k, gUnm1, gVnm1, myThid )               CALL OBCS_APPLY_UV( bi, bj, k, gU, gV, myThid )
406           END IF             ENDIF
407  #endif   /* ALLOW_OBCS */  #endif   /* ALLOW_OBCS */
408    
 #ifdef   ALLOW_AUTODIFF_TAMC  
 #ifdef   INCLUDE_CD_CODE  
          ELSE  
            DO j=1-OLy,sNy+OLy  
              DO i=1-OLx,sNx+OLx  
                guCD(i,j,k,bi,bj) = 0.0  
                gvCD(i,j,k,bi,bj) = 0.0  
              END DO  
            END DO  
 #endif   /* INCLUDE_CD_CODE */  
 #endif   /* ALLOW_AUTODIFF_TAMC */  
409           ENDIF           ENDIF
410    
411    
412  C--     end of dynamics k loop (1:Nr)  C--     end of dynamics k loop (1:Nr)
413          ENDDO          ENDDO
414    
   
   
415  C--     Implicit viscosity  C--     Implicit viscosity
416          IF (implicitViscosity.AND.momStepping) THEN          IF (implicitViscosity.AND.momStepping) THEN
417  #ifdef    ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
418            idkey = iikey + 3  CADJ STORE KappaRU(:,:,:) = comlev1_bibj , key=idynkey, byte=isbyte
419  CADJ STORE gUNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE gU(:,:,:,bi,bj) = comlev1_bibj , key=idynkey, byte=isbyte
420  #endif    /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
421            CALL IMPLDIFF(            CALL IMPLDIFF(
422       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
423       I         deltaTmom, KappaRU,recip_HFacW,       I         deltaTmom, KappaRU,recip_HFacW,
424       U         gUNm1,       U         gU,
425       I         myThid )       I         myThid )
426  #ifdef    ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
427            idkey = iikey + 4  CADJ STORE KappaRV(:,:,:) = comlev1_bibj , key=idynkey, byte=isbyte
428  CADJ STORE gVNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE gV(:,:,:,bi,bj) = comlev1_bibj , key=idynkey, byte=isbyte
429  #endif    /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
430            CALL IMPLDIFF(            CALL IMPLDIFF(
431       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
432       I         deltaTmom, KappaRV,recip_HFacS,       I         deltaTmom, KappaRV,recip_HFacS,
433       U         gVNm1,       U         gV,
434       I         myThid )       I         myThid )
435    
436  #ifdef   ALLOW_OBCS  #ifdef   ALLOW_OBCS
437  C--      Apply open boundary conditions  C--      Apply open boundary conditions
438           IF (useOBCS) THEN           IF (useOBCS) THEN
439             DO K=1,Nr             DO K=1,Nr
440               CALL OBCS_APPLY_UV( bi, bj, k, gUnm1, gVnm1, myThid )               CALL OBCS_APPLY_UV( bi, bj, k, gU, gV, myThid )
441             ENDDO             ENDDO
442           END IF           END IF
443  #endif   /* ALLOW_OBCS */  #endif   /* ALLOW_OBCS */
444    
445  #ifdef    INCLUDE_CD_CODE  #ifdef    ALLOW_CD_CODE
446  #ifdef    ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
447            idkey = iikey + 5  CADJ STORE vVelD(:,:,:,bi,bj) = comlev1_bibj , key=idynkey, byte=isbyte
 CADJ STORE vVelD(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  
448  #endif    /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
449            CALL IMPLDIFF(            CALL IMPLDIFF(
450       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
# Line 717  CADJ STORE vVelD(:,:,:,bi,bj) = comlev1_ Line 452  CADJ STORE vVelD(:,:,:,bi,bj) = comlev1_
452       U         vVelD,       U         vVelD,
453       I         myThid )       I         myThid )
454  #ifdef    ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
455            idkey = iikey + 6  CADJ STORE uVelD(:,:,:,bi,bj) = comlev1_bibj , key=idynkey, byte=isbyte
 CADJ STORE uVelD(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  
456  #endif    /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
457            CALL IMPLDIFF(            CALL IMPLDIFF(
458       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
459       I         deltaTmom, KappaRV,recip_HFacS,       I         deltaTmom, KappaRV,recip_HFacS,
460       U         uVelD,       U         uVelD,
461       I         myThid )       I         myThid )
462  #endif    /* INCLUDE_CD_CODE */  #endif    /* ALLOW_CD_CODE */
463  C--     End If implicitViscosity.AND.momStepping  C--     End If implicitViscosity.AND.momStepping
464          ENDIF          ENDIF
465    
 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_CUMUL_1T(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 */  
   
466         ENDDO         ENDDO
467        ENDDO        ENDDO
468    
469    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 ( debugLevel .GE. debLevB ) THEN
479           CALL DEBUG_STATS_RL(1,EtaN,'EtaN (DYNAMICS)',myThid)
480           CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (DYNAMICS)',myThid)
481           CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (DYNAMICS)',myThid)
482           CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (DYNAMICS)',myThid)
483           CALL DEBUG_STATS_RL(Nr,theta,'Theta (DYNAMICS)',myThid)
484           CALL DEBUG_STATS_RL(Nr,salt,'Salt (DYNAMICS)',myThid)
485           CALL DEBUG_STATS_RL(Nr,Gu,'Gu (DYNAMICS)',myThid)
486           CALL DEBUG_STATS_RL(Nr,Gv,'Gv (DYNAMICS)',myThid)
487           CALL DEBUG_STATS_RL(Nr,Gt,'Gt (DYNAMICS)',myThid)
488           CALL DEBUG_STATS_RL(Nr,Gs,'Gs (DYNAMICS)',myThid)
489           CALL DEBUG_STATS_RL(Nr,GuNm1,'GuNm1 (DYNAMICS)',myThid)
490           CALL DEBUG_STATS_RL(Nr,GvNm1,'GvNm1 (DYNAMICS)',myThid)
491           CALL DEBUG_STATS_RL(Nr,GtNm1,'GtNm1 (DYNAMICS)',myThid)
492           CALL DEBUG_STATS_RL(Nr,GsNm1,'GsNm1 (DYNAMICS)',myThid)
493          ENDIF
494    #endif
495    
496        RETURN        RETURN
497        END        END
Legend:
Removed from v.1.68  
changed lines
  Added in v.1.103
  ViewVC Help
Powered by ViewVC 1.1.22