/[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.69 by adcroft, Wed Jun 6 14:55:45 2001 UTC revision 1.133 by heimbach, Wed May 31 19:53:28 2006 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    #ifdef ALLOW_OBCS
7    # include "OBCS_OPTIONS.h"
8    #endif
9    
10    #undef DYNAMICS_GUGV_EXCH_CHECK
11    
12    CBOP
13    C     !ROUTINE: DYNAMICS
14    C     !INTERFACE:
15        SUBROUTINE DYNAMICS(myTime, myIter, myThid)        SUBROUTINE DYNAMICS(myTime, myIter, myThid)
16  C     /==========================================================\  C     !DESCRIPTION: \bv
17  C     | SUBROUTINE DYNAMICS                                      |  C     *==========================================================*
18  C     | o Controlling routine for the explicit part of the model |  C     | SUBROUTINE DYNAMICS                                      
19  C     |   dynamics.                                              |  C     | o Controlling routine for the explicit part of the model  
20  C     |==========================================================|  C     |   dynamics.                                              
21  C     | This routine evaluates the "dynamics" terms for each     |  C     *==========================================================*
22  C     | block of ocean in turn. Because the blocks of ocean have |  C     | This routine evaluates the "dynamics" terms for each      
23  C     | overlap regions they are independent of one another.     |  C     | block of ocean in turn. Because the blocks of ocean have  
24  C     | If terms involving lateral integrals are needed in this  |  C     | overlap regions they are independent of one another.      
25  C     | routine care will be needed. Similarly finite-difference |  C     | If terms involving lateral integrals are needed in this  
26  C     | operations with stencils wider than the overlap region   |  C     | routine care will be needed. Similarly finite-difference  
27  C     | require special consideration.                           |  C     | operations with stencils wider than the overlap region    
28  C     | Notes                                                    |  C     | require special consideration.                            
29  C     | =====                                                    |  C     | The algorithm...
30  C     | C*P* comments indicating place holders for which code is |  C     |
31  C     |      presently being developed.                          |  C     | "Correction Step"
32  C     \==========================================================/  C     | =================
33    C     | Here we update the horizontal velocities with the surface
34    C     | pressure such that the resulting flow is either consistent
35    C     | with the free-surface evolution or the rigid-lid:
36    C     |   U[n] = U* + dt x d/dx P
37    C     |   V[n] = V* + dt x d/dy P
38    C     |   W[n] = W* + dt x d/dz P  (NH mode)
39    C     |
40    C     | "Calculation of Gs"
41    C     | ===================
42    C     | This is where all the accelerations and tendencies (ie.
43    C     | physics, parameterizations etc...) are calculated
44    C     |   rho = rho ( theta[n], salt[n] )
45    C     |   b   = b(rho, theta)
46    C     |   K31 = K31 ( rho )
47    C     |   Gu[n] = Gu( u[n], v[n], wVel, b, ... )
48    C     |   Gv[n] = Gv( u[n], v[n], wVel, b, ... )
49    C     |   Gt[n] = Gt( theta[n], u[n], v[n], wVel, K31, ... )
50    C     |   Gs[n] = Gs( salt[n], u[n], v[n], wVel, K31, ... )
51    C     |
52    C     | "Time-stepping" or "Prediction"
53    C     | ================================
54    C     | The models variables are stepped forward with the appropriate
55    C     | time-stepping scheme (currently we use Adams-Bashforth II)
56    C     | - For momentum, the result is always *only* a "prediction"
57    C     | in that the flow may be divergent and will be "corrected"
58    C     | later with a surface pressure gradient.
59    C     | - Normally for tracers the result is the new field at time
60    C     | level [n+1} *BUT* in the case of implicit diffusion the result
61    C     | is also *only* a prediction.
62    C     | - We denote "predictors" with an asterisk (*).
63    C     |   U* = U[n] + dt x ( 3/2 Gu[n] - 1/2 Gu[n-1] )
64    C     |   V* = V[n] + dt x ( 3/2 Gv[n] - 1/2 Gv[n-1] )
65    C     |   theta[n+1] = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
66    C     |   salt[n+1] = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
67    C     | With implicit diffusion:
68    C     |   theta* = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
69    C     |   salt* = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
70    C     |   (1 + dt * K * d_zz) theta[n] = theta*
71    C     |   (1 + dt * K * d_zz) salt[n] = salt*
72    C     |
73    C     *==========================================================*
74    C     \ev
75    C     !USES:
76        IMPLICIT NONE        IMPLICIT NONE
   
77  C     == Global variables ===  C     == Global variables ===
78  #include "SIZE.h"  #include "SIZE.h"
79  #include "EEPARAMS.h"  #include "EEPARAMS.h"
80  #include "PARAMS.h"  #include "PARAMS.h"
81  #include "DYNVARS.h"  #include "DYNVARS.h"
82    #ifdef ALLOW_CD_CODE
83    #include "CD_CODE_VARS.h"
84    #endif
85  #include "GRID.h"  #include "GRID.h"
   
86  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
87  # include "tamc.h"  # include "tamc.h"
88  # include "tamc_keys.h"  # include "tamc_keys.h"
89  # include "FFIELDS.h"  # include "FFIELDS.h"
90    # include "EOS.h"
91  # ifdef ALLOW_KPP  # ifdef ALLOW_KPP
92  #  include "KPP.h"  #  include "KPP.h"
93  # endif  # endif
94  # ifdef ALLOW_GMREDI  # ifdef ALLOW_PTRACERS
95  #  include "GMREDI.h"  #  include "PTRACERS_SIZE.h"
96    #  include "PTRACERS.h"
97    # endif
98    # ifdef ALLOW_OBCS
99    #  include "OBCS.h"
100    #  ifdef ALLOW_PTRACERS
101    #   include "OBCS_PTRACERS.h"
102    #  endif
103    # endif
104    # ifdef ALLOW_MOM_FLUXFORM
105    #  include "MOM_FLUXFORM.h"
106  # endif  # endif
107  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
108    
109  #ifdef ALLOW_TIMEAVE  C     !CALLING SEQUENCE:
110  #include "TIMEAVE_STATV.h"  C     DYNAMICS()
111  #endif  C      |
112    C      |-- CALC_EP_FORCING
113    C      |
114    C      |-- CALC_GRAD_PHI_SURF
115    C      |
116    C      |-- CALC_VISCOSITY
117    C      |
118    C      |-- CALC_PHI_HYD  
119    C      |
120    C      |-- MOM_FLUXFORM  
121    C      |
122    C      |-- MOM_VECINV    
123    C      |
124    C      |-- TIMESTEP      
125    C      |
126    C      |-- OBCS_APPLY_UV
127    C      |
128    C      |-- MOM_U_IMPLICIT_R      
129    C      |-- MOM_V_IMPLICIT_R      
130    C      |
131    C      |-- IMPLDIFF      
132    C      |
133    C      |-- OBCS_APPLY_UV
134    C      |
135    C      |-- CALC_GW
136    C      |
137    C      |-- DIAGNOSTICS_FILL
138    C      |-- DEBUG_STATS_RL
139    
140    C     !INPUT/OUTPUT PARAMETERS:
141  C     == Routine arguments ==  C     == Routine arguments ==
142  C     myTime - Current time in simulation  C     myTime - Current time in simulation
143  C     myIter - Current iteration number in simulation  C     myIter - Current iteration number in simulation
# Line 54  C     myThid - Thread number for this in Line 146  C     myThid - Thread number for this in
146        INTEGER myIter        INTEGER myIter
147        INTEGER myThid        INTEGER myThid
148    
149    C     !LOCAL VARIABLES:
150  C     == Local variables  C     == Local variables
151  C     xA, yA                 - Per block temporaries holding face areas  C     fVer[UV]               o fVer: Vertical flux term - note fVer
152  C     uTrans, vTrans, rTrans - Per block temporaries holding flow  C                                    is "pipelined" in the vertical
153  C                              transport  C                                    so we need an fVer for each
154  C                              o uTrans: Zonal transport  C                                    variable.
155  C                              o vTrans: Meridional transport  C     phiHydC    :: hydrostatic potential anomaly at cell center
156  C                              o rTrans: Vertical transport  C                   In z coords phiHyd is the hydrostatic potential
157  C     maskUp                   o maskUp: land/water mask for W points  C                      (=pressure/rho0) anomaly
158  C     fVer[STUV]               o fVer: Vertical flux term - note fVer  C                   In p coords phiHyd is the geopotential height anomaly.
159  C                                      is "pipelined" in the vertical  C     phiHydF    :: hydrostatic potential anomaly at middle between 2 centers
160  C                                      so we need an fVer for each  C     dPhiHydX,Y :: Gradient (X & Y directions) of hydrostatic potential anom.
161  C                                      variable.  C     phiSurfX,  ::  gradient of Surface potential (Pressure/rho, ocean)
162  C     rhoK, rhoKM1   - Density at current level, and level above  C     phiSurfY             or geopotential (atmos) in X and Y direction
163  C     phiHyd         - Hydrostatic part of the potential phiHydi.  C     guDissip   :: dissipation tendency (all explicit terms), u component
164  C                      In z coords phiHydiHyd is the hydrostatic  C     gvDissip   :: dissipation tendency (all explicit terms), v component
 C                      Potential (=pressure/rho0) anomaly  
 C                      In p coords phiHydiHyd is the geopotential  
 C                      surface height anomaly.  
 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).  
165  C     iMin, iMax     - Ranges and sub-block indices on which calculations  C     iMin, iMax     - Ranges and sub-block indices on which calculations
166  C     jMin, jMax       are applied.  C     jMin, jMax       are applied.
167  C     bi, bj  C     bi, bj
168  C     k, kup,        - Index for layer above and below. kup and kDown  C     k, kup,        - Index for layer above and below. kup and kDown
169  C     kDown, km1       are switched with layer to be the appropriate  C     kDown, km1       are switched with layer to be the appropriate
170  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)  
171        _RL fVerU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
172        _RL fVerV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
173        _RL phiHyd  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL phiHydF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
174        _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiHydC (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
175        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL dPhiHydX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
176          _RL dPhiHydY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
177        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
178        _RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
179        _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL guDissip(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
180        _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL gvDissip(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
181        _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
182        _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)  
183    
184        INTEGER iMin, iMax        INTEGER iMin, iMax
185        INTEGER jMin, jMax        INTEGER jMin, jMax
186        INTEGER bi, bj        INTEGER bi, bj
187        INTEGER i, j        INTEGER i, j
188        INTEGER k, km1, kup, kDown        INTEGER k, km1, kp1, kup, kDown
189    
190    #ifdef ALLOW_DIAGNOSTICS
191          _RL tmpFac
192    #endif /* ALLOW_DIAGNOSTICS */
193    
 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)  
194    
195  C---    The algorithm...  C---    The algorithm...
196  C  C
# Line 165  C         salt* = salt[n] + dt x ( 3/2 G Line 235  C         salt* = salt[n] + dt x ( 3/2 G
235  C         (1 + dt * K * d_zz) theta[n] = theta*  C         (1 + dt * K * d_zz) theta[n] = theta*
236  C         (1 + dt * K * d_zz) salt[n] = salt*  C         (1 + dt * K * d_zz) salt[n] = salt*
237  C---  C---
238    CEOP
239    
240  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_DEBUG
241  C--   dummy statement to end declaration part        IF ( debugLevel .GE. debLevB )
242        ikey = 1       &   CALL DEBUG_ENTER( 'DYNAMICS', myThid )
243  #endif /* ALLOW_AUTODIFF_TAMC */  #endif
   
 C--   Set up work arrays with valid (i.e. not NaN) values  
 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  
244    
245    C-- Call to routine for calculation of
246    C   Eliassen-Palm-flux-forced U-tendency,
247    C   if desired:
248    #ifdef INCLUDE_EP_FORCING_CODE
249          CALL CALC_EP_FORCING(myThid)
250    #endif
251    
252  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
253  C--   HPF directive to help TAMC  C--   HPF directive to help TAMC
# Line 207  CHPF$ INDEPENDENT Line 258  CHPF$ INDEPENDENT
258    
259  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
260  C--    HPF directive to help TAMC  C--    HPF directive to help TAMC
261  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS,fVerU,fVerV  CHPF$  INDEPENDENT, NEW (fVerU,fVerV
262  CHPF$&                  ,phiHyd,utrans,vtrans,xA,yA  CHPF$&                  ,phiHydF
263  CHPF$&                  ,KappaRT,KappaRS,KappaRU,KappaRV  CHPF$&                  ,KappaRU,KappaRV
264  CHPF$&                  )  CHPF$&                  )
265  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
266    
# Line 218  CHPF$&                  ) Line 269  CHPF$&                  )
269  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
270            act1 = bi - myBxLo(myThid)            act1 = bi - myBxLo(myThid)
271            max1 = myBxHi(myThid) - myBxLo(myThid) + 1            max1 = myBxHi(myThid) - myBxLo(myThid) + 1
   
272            act2 = bj - myByLo(myThid)            act2 = bj - myByLo(myThid)
273            max2 = myByHi(myThid) - myByLo(myThid) + 1            max2 = myByHi(myThid) - myByLo(myThid) + 1
   
274            act3 = myThid - 1            act3 = myThid - 1
275            max3 = nTx*nTy            max3 = nTx*nTy
   
276            act4 = ikey_dynamics - 1            act4 = ikey_dynamics - 1
277              idynkey = (act1 + 1) + act2*max1
           ikey = (act1 + 1) + act2*max1  
278       &                      + act3*max1*max2       &                      + act3*max1*max2
279       &                      + act4*max1*max2*max3       &                      + act4*max1*max2*max3
280  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
281    
282  C--     Set up work arrays that need valid initial values  C--   Set up work arrays with valid (i.e. not NaN) values
283          DO j=1-OLy,sNy+OLy  C     These inital values do not alter the numerical results. They
284           DO i=1-OLx,sNx+OLx  C     just ensure that all memory references are to valid floating
285            rTrans(i,j)   = 0. _d 0  C     point numbers. This prevents spurious hardware signals due to
286            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  
287    
288          DO k=1,Nr          DO k=1,Nr
289           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
290            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
291  C This is currently also used by IVDC and Diagnostics             KappaRU(i,j,k) = 0. _d 0
292             ConvectCount(i,j,k) = 0.             KappaRV(i,j,k) = 0. _d 0
293             KappaRT(i,j,k) = 0. _d 0  #ifdef ALLOW_AUTODIFF_TAMC
294             KappaRS(i,j,k) = 0. _d 0  cph(
295    c--   need some re-initialisation here to break dependencies
296    cph)
297               gU(i,j,k,bi,bj) = 0. _d 0
298               gV(i,j,k,bi,bj) = 0. _d 0
299    #endif
300            ENDDO            ENDDO
301           ENDDO           ENDDO
302          ENDDO          ENDDO
303            DO j=1-OLy,sNy+OLy
304          iMin = 1-OLx+1           DO i=1-OLx,sNx+OLx
305          iMax = sNx+OLx            fVerU  (i,j,1) = 0. _d 0
306          jMin = 1-OLy+1            fVerU  (i,j,2) = 0. _d 0
307          jMax = sNy+OLy            fVerV  (i,j,1) = 0. _d 0
308              fVerV  (i,j,2) = 0. _d 0
309              phiHydF (i,j)  = 0. _d 0
310  #ifdef ALLOW_AUTODIFF_TAMC            phiHydC (i,j)  = 0. _d 0
311  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte            dPhiHydX(i,j)  = 0. _d 0
312  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte            dPhiHydY(i,j)  = 0. _d 0
313  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte            phiSurfX(i,j)  = 0. _d 0
314  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte            phiSurfY(i,j)  = 0. _d 0
315  #endif /* ALLOW_AUTODIFF_TAMC */            guDissip(i,j)  = 0. _d 0
316              gvDissip(i,j)  = 0. _d 0
317  C--     Start of diagnostic loop  #ifdef ALLOW_AUTODIFF_TAMC
318          DO k=Nr,1,-1  cph(
319    c--   need some re-initialisation here to break dependencies
320  #ifdef ALLOW_AUTODIFF_TAMC  cph)
321  C? Patrick, is this formula correct now that we change the loop range?  # ifdef NONLIN_FRSURF
322  C? Do we still need this?  #  ifndef DISABLE_RSTAR_CODE
323  cph kkey formula corrected.            dWtransC(i,j,bi,bj)  = 0. _d 0
324  cph Needed for rhok, rhokm1, in the case useGMREDI.            dWtransU(i,j,bi,bj)  = 0. _d 0
325           kkey = (ikey-1)*Nr + k            dWtransV(i,j,bi,bj)  = 0. _d 0
326  CADJ STORE rhokm1(:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte  #  endif
327  CADJ STORE rhok  (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte  # endif /* NONLIN_FRSURF */
 #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  
328  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
329               CALL FIND_RHO(           ENDDO
      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)  
330          ENDDO          ENDDO
331    
332  #ifdef ALLOW_AUTODIFF_TAMC  C--     Start computation of dynamics
333  cph avoids recomputation of integrate_for_w          iMin = 0
334  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte          iMax = sNx+1
335  #endif /* ALLOW_AUTODIFF_TAMC */          jMin = 0
336            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  
337    
338  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
339  CADJ STORE Kwx(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE wvel (:,:,:,bi,bj) =
340  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  
341  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
342    
343  #endif  /* ALLOW_GMREDI */  C--     Explicit part of the Surface Potentiel Gradient (add in TIMESTEP)
344    C       (note: this loop will be replaced by CALL CALC_GRAD_ETA)
345  #ifdef  ALLOW_KPP          IF (implicSurfPress.NE.1.) THEN
346  C--     Compute KPP mixing coefficients            CALL CALC_GRAD_PHI_SURF(
347          IF (useKPP) THEN       I         bi,bj,iMin,iMax,jMin,jMax,
348            CALL KPP_CALC(       I         etaN,
349       I                  bi, bj, myTime, myThid )       O         phiSurfX,phiSurfY,
350  #ifdef ALLOW_AUTODIFF_TAMC       I         myThid )                        
         ELSE  
           CALL KPP_CALC_DUMMY(  
      I                  bi, bj, myTime, myThid )  
 #endif /* ALLOW_AUTODIFF_TAMC */  
351          ENDIF          ENDIF
352    
353  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
354  CADJ STORE KPPghat   (:,:,:,bi,bj)  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=idynkey, byte=isbyte
355  CADJ &   , KPPviscAz (:,:,:,bi,bj)  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=idynkey, byte=isbyte
356  CADJ &   , KPPdiffKzT(:,:,:,bi,bj)  #ifdef ALLOW_KPP
357  CADJ &   , KPPdiffKzS(:,:,:,bi,bj)  CADJ STORE KPPviscAz (:,:,:,bi,bj)
358  CADJ &   , KPPfrac   (:,:  ,bi,bj)  CADJ &                 = comlev1_bibj, key=idynkey, byte=isbyte
359  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  
360  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
361    
362  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL
363  C--      Calculate the total vertical diffusivity  C--      Calculate the total vertical diffusivity
364           CALL CALC_DIFFUSIVITY(          DO k=1,Nr
365             CALL CALC_VISCOSITY(
366       I        bi,bj,iMin,iMax,jMin,jMax,k,       I        bi,bj,iMin,iMax,jMin,jMax,k,
367       I        maskUp,       O        KappaRU,KappaRV,
      O        KappaRT,KappaRS,KappaRU,KappaRV,  
368       I        myThid)       I        myThid)
369           ENDDO
370  #endif  #endif
371    
 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  
372  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
373  CADJ STORE gTNm1(:,:,k,bi,bj) = comlev1_bibj_k  CADJ STORE KappaRU(:,:,:)
374  CADJ &   , key = kkey, byte = isbyte  CADJ &     = comlev1_bibj, key=idynkey, byte=isbyte
375    CADJ STORE KappaRV(:,:,:)
376    CADJ &     = comlev1_bibj, key=idynkey, byte=isbyte
377  #endif /* ALLOW_AUTODIFF_TAMC */  #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_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  
378    
379  C--     Start of dynamics loop  C--     Start of dynamics loop
380          DO k=1,Nr          DO k=1,Nr
# Line 617  C--       kup    Cycles through 1,2 to p Line 384  C--       kup    Cycles through 1,2 to p
384  C--       kDown  Cycles through 2,1 to point to current layer  C--       kDown  Cycles through 2,1 to point to current layer
385    
386            km1  = MAX(1,k-1)            km1  = MAX(1,k-1)
387              kp1  = MIN(k+1,Nr)
388            kup  = 1+MOD(k+1,2)            kup  = 1+MOD(k+1,2)
389            kDown= 1+MOD(k,2)            kDown= 1+MOD(k,2)
390    
391    #ifdef ALLOW_AUTODIFF_TAMC
392             kkey = (idynkey-1)*Nr + k
393    c
394    CADJ STORE totphihyd (:,:,k,bi,bj)
395    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
396    CADJ STORE theta (:,:,k,bi,bj)
397    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
398    CADJ STORE salt  (:,:,k,bi,bj)
399    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
400    CADJ STORE gt(:,:,k,bi,bj)
401    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
402    CADJ STORE gs(:,:,k,bi,bj)
403    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
404    # ifdef NONLIN_FRSURF
405    cph-test
406    CADJ STORE  phiHydC (:,:)
407    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
408    CADJ STORE  phiHydF (:,:)
409    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
410    CADJ STORE  gudissip (:,:)
411    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
412    CADJ STORE  gvdissip (:,:)
413    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
414    CADJ STORE  fVerU (:,:,:)
415    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
416    CADJ STORE  fVerV (:,:,:)
417    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
418    CADJ STORE gu(:,:,k,bi,bj)
419    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
420    CADJ STORE gv(:,:,k,bi,bj)
421    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
422    CADJ STORE gunm1(:,:,k,bi,bj)
423    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
424    CADJ STORE gvnm1(:,:,k,bi,bj)
425    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
426    #   ifndef DISABLE_RSTAR_CODE
427    CADJ STORE dwtransc(:,:,bi,bj)
428    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
429    CADJ STORE dwtransu(:,:,bi,bj)
430    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
431    CADJ STORE dwtransv(:,:,bi,bj)
432    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
433    #   endif
434    #  ifdef ALLOW_CD_CODE
435    CADJ STORE unm1(:,:,k,bi,bj)
436    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
437    CADJ STORE vnm1(:,:,k,bi,bj)
438    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
439    CADJ STORE uVelD(:,:,k,bi,bj)
440    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
441    CADJ STORE vVelD(:,:,k,bi,bj)
442    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
443    #  endif
444    # endif
445    #endif /* ALLOW_AUTODIFF_TAMC */
446    
447  C--      Integrate hydrostatic balance for phiHyd with BC of  C--      Integrate hydrostatic balance for phiHyd with BC of
448  C        phiHyd(z=0)=0  C        phiHyd(z=0)=0
449  C        distinguishe between Stagger and Non Stagger time stepping           IF ( implicitIntGravWave ) THEN
          IF (staggerTimeStep) THEN  
450             CALL CALC_PHI_HYD(             CALL CALC_PHI_HYD(
451       I        bi,bj,iMin,iMax,jMin,jMax,k,       I        bi,bj,iMin,iMax,jMin,jMax,k,
452       I        gTnm1, gSnm1,       I        gT, gS,
453       U        phiHyd,       U        phiHydF,
454       I        myThid )       O        phiHydC, dPhiHydX, dPhiHydY,
455         I        myTime, myIter, myThid )
456           ELSE           ELSE
457             CALL CALC_PHI_HYD(             CALL CALC_PHI_HYD(
458       I        bi,bj,iMin,iMax,jMin,jMax,k,       I        bi,bj,iMin,iMax,jMin,jMax,k,
459       I        theta, salt,       I        theta, salt,
460       U        phiHyd,       U        phiHydF,
461       I        myThid )       O        phiHydC, dPhiHydX, dPhiHydY,
462         I        myTime, myIter, myThid )
463           ENDIF           ENDIF
464    
465  C--      Calculate accelerations in the momentum equations (gU, gV, ...)  C--      Calculate accelerations in the momentum equations (gU, gV, ...)
466  C        and step forward storing the result in gUnm1, gVnm1, etc...  C        and step forward storing the result in gU, gV, etc...
467           IF ( momStepping ) THEN           IF ( momStepping ) THEN
468             CALL CALC_MOM_RHS(             IF (.NOT. vectorInvariantMomentum) THEN
469    #ifdef ALLOW_MOM_FLUXFORM
470    C
471    # ifdef ALLOW_AUTODIFF_TAMC
472    #  ifdef NONLIN_FRSURF
473    #   ifndef DISABLE_RSTAR_CODE
474    CADJ STORE dwtransc(:,:,bi,bj)
475    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
476    CADJ STORE dwtransu(:,:,bi,bj)
477    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
478    CADJ STORE dwtransv(:,:,bi,bj)
479    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
480    #   endif
481    #  endif
482    # endif /* ALLOW_AUTODIFF_TAMC */
483    C
484                  CALL MOM_FLUXFORM(
485       I         bi,bj,iMin,iMax,jMin,jMax,k,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax,k,kup,kDown,
486       I         phiHyd,KappaRU,KappaRV,       I         KappaRU, KappaRV,
487       U         fVerU, fVerV,       U         fVerU, fVerV,
488       I         myTime, myThid)       O         guDissip, gvDissip,
489         I         myTime, myIter, myThid)
490    #endif
491               ELSE
492    #ifdef ALLOW_MOM_VECINV
493    C
494    # ifdef ALLOW_AUTODIFF_TAMC
495    #  ifdef NONLIN_FRSURF
496    CADJ STORE fVerU(:,:,:)
497    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
498    CADJ STORE fVerV(:,:,:)
499    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
500    #  endif
501    # endif /* ALLOW_AUTODIFF_TAMC */
502    C
503                 CALL MOM_VECINV(
504         I         bi,bj,iMin,iMax,jMin,jMax,k,kup,kDown,
505         I         KappaRU, KappaRV,
506         U         fVerU, fVerV,
507         O         guDissip, gvDissip,
508         I         myTime, myIter, myThid)
509    #endif
510               ENDIF
511    C
512             CALL TIMESTEP(             CALL TIMESTEP(
513       I         bi,bj,iMin,iMax,jMin,jMax,k,       I         bi,bj,iMin,iMax,jMin,jMax,k,
514       I         phiHyd, phiSurfX, phiSurfY,       I         dPhiHydX,dPhiHydY, phiSurfX, phiSurfY,
515       I         myIter, myThid)       I         guDissip, gvDissip,
516         I         myTime, myIter, myThid)
517    
518  #ifdef   ALLOW_OBCS  #ifdef   ALLOW_OBCS
519  C--      Apply open boundary conditions  C--      Apply open boundary conditions
520           IF (useOBCS) THEN             IF (useOBCS) THEN
521             CALL OBCS_APPLY_UV( bi, bj, k, gUnm1, gVnm1, myThid )               CALL OBCS_APPLY_UV( bi, bj, k, gU, gV, myThid )
522           END IF             ENDIF
523  #endif   /* ALLOW_OBCS */  #endif   /* ALLOW_OBCS */
524    
 #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 */  
525           ENDIF           ENDIF
526    
527    
528  C--     end of dynamics k loop (1:Nr)  C--     end of dynamics k loop (1:Nr)
529          ENDDO          ENDDO
530    
531    C--     Implicit Vertical advection & viscosity
532    #if (defined (INCLUDE_IMPLVERTADV_CODE) && defined (ALLOW_MOM_COMMON))
533  C--     Implicit viscosity          IF ( momImplVertAdv ) THEN
534          IF (implicitViscosity.AND.momStepping) THEN            CALL MOM_U_IMPLICIT_R( kappaRU,
535         I                           bi, bj, myTime, myIter, myThid )
536              CALL MOM_V_IMPLICIT_R( kappaRV,
537         I                           bi, bj, myTime, myIter, myThid )
538            ELSEIF ( implicitViscosity ) THEN
539    #else /* INCLUDE_IMPLVERTADV_CODE */
540            IF     ( implicitViscosity ) THEN
541    #endif /* INCLUDE_IMPLVERTADV_CODE */
542  #ifdef    ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
543            idkey = iikey + 3  CADJ STORE KappaRU(:,:,:) = comlev1_bibj , key=idynkey, byte=isbyte
544  CADJ STORE gUNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE gU(:,:,:,bi,bj) = comlev1_bibj , key=idynkey, byte=isbyte
545  #endif    /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
546            CALL IMPLDIFF(            CALL IMPLDIFF(
547       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
548       I         deltaTmom, KappaRU,recip_HFacW,       I         -1, KappaRU,recip_HFacW,
549       U         gUNm1,       U         gU,
550       I         myThid )       I         myThid )
551  #ifdef    ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
552            idkey = iikey + 4  CADJ STORE KappaRV(:,:,:) = comlev1_bibj , key=idynkey, byte=isbyte
553  CADJ STORE gVNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE gV(:,:,:,bi,bj) = comlev1_bibj , key=idynkey, byte=isbyte
554  #endif    /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
555            CALL IMPLDIFF(            CALL IMPLDIFF(
556       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
557       I         deltaTmom, KappaRV,recip_HFacS,       I         -2, KappaRV,recip_HFacS,
558       U         gVNm1,       U         gV,
559       I         myThid )       I         myThid )
560            ENDIF
561    
562  #ifdef   ALLOW_OBCS  #ifdef   ALLOW_OBCS
563  C--      Apply open boundary conditions  C--      Apply open boundary conditions
564           IF (useOBCS) THEN          IF ( useOBCS .AND.(implicitViscosity.OR.momImplVertAdv) ) THEN
565             DO K=1,Nr             DO K=1,Nr
566               CALL OBCS_APPLY_UV( bi, bj, k, gUnm1, gVnm1, myThid )               CALL OBCS_APPLY_UV( bi, bj, k, gU, gV, myThid )
567             ENDDO             ENDDO
568           END IF          ENDIF
569  #endif   /* ALLOW_OBCS */  #endif   /* ALLOW_OBCS */
570    
571  #ifdef    INCLUDE_CD_CODE  #ifdef    ALLOW_CD_CODE
572            IF (implicitViscosity.AND.useCDscheme) THEN
573  #ifdef    ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
574            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  
575  #endif    /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
576            CALL IMPLDIFF(            CALL IMPLDIFF(
577       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
578       I         deltaTmom, KappaRU,recip_HFacW,       I         0, KappaRU,recip_HFacW,
579       U         vVelD,       U         vVelD,
580       I         myThid )       I         myThid )
581  #ifdef    ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
582            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  
583  #endif    /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
584            CALL IMPLDIFF(            CALL IMPLDIFF(
585       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
586       I         deltaTmom, KappaRV,recip_HFacS,       I         0, KappaRV,recip_HFacS,
587       U         uVelD,       U         uVelD,
588       I         myThid )       I         myThid )
 #endif    /* INCLUDE_CD_CODE */  
 C--     End If implicitViscosity.AND.momStepping  
589          ENDIF          ENDIF
590    #endif    /* ALLOW_CD_CODE */
591    C--     End implicit Vertical advection & viscosity
592    
 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 */  
   
593         ENDDO         ENDDO
594        ENDDO        ENDDO
595    
596  #ifndef EXCLUDE_DEBUGMODE  #ifdef ALLOW_OBCS
597          IF (useOBCS) THEN
598           CALL OBCS_PRESCRIBE_EXCHANGES(myThid)
599          ENDIF
600    #endif
601    
602    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
603    
604    #ifdef ALLOW_NONHYDROSTATIC
605    C--   Step forward W field in N-H algorithm
606          IF ( nonHydrostatic ) THEN
607    #ifdef ALLOW_DEBUG
608             IF ( debugLevel .GE. debLevB )
609         &     CALL DEBUG_CALL('CALC_GW', myThid )
610    #endif
611             CALL TIMER_START('CALC_GW          [DYNAMICS]',myThid)
612             CALL CALC_GW( myTime, myIter, myThid )
613          ENDIF
614          IF ( nonHydrostatic.OR.implicitIntGravWave )
615         &   CALL TIMESTEP_WVEL( myTime, myIter, myThid )
616          IF ( nonHydrostatic )
617         &   CALL TIMER_STOP ('CALC_GW          [DYNAMICS]',myThid)
618    #endif
619    
620    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
621    
622    Cml(
623    C     In order to compare the variance of phiHydLow of a p/z-coordinate
624    C     run with etaH of a z/p-coordinate run the drift of phiHydLow
625    C     has to be removed by something like the following subroutine:
626    C      CALL REMOVE_MEAN_RL( 1, phiHydLow, maskH, maskH, rA, drF,
627    C     &                'phiHydLow', myThid )
628    Cml)
629    
630    #ifdef ALLOW_DIAGNOSTICS
631          IF ( useDiagnostics ) THEN
632    
633           CALL DIAGNOSTICS_FILL(totPhihyd,'PHIHYD  ',0,Nr,0,1,1,myThid)
634           CALL DIAGNOSTICS_FILL(phiHydLow,'PHIBOT  ',0, 1,0,1,1,myThid)
635    
636           tmpFac = 1. _d 0
637           CALL DIAGNOSTICS_SCALE_FILL(totPhihyd,tmpFac,2,
638         &                                 'PHIHYDSQ',0,Nr,0,1,1,myThid)
639    
640           CALL DIAGNOSTICS_SCALE_FILL(phiHydLow,tmpFac,2,
641         &                                 'PHIBOTSQ',0, 1,0,1,1,myThid)
642    
643          ENDIF
644    #endif /* ALLOW_DIAGNOSTICS */
645          
646    #ifdef ALLOW_DEBUG
647          If ( debugLevel .GE. debLevB ) THEN
648         CALL DEBUG_STATS_RL(1,EtaN,'EtaN (DYNAMICS)',myThid)         CALL DEBUG_STATS_RL(1,EtaN,'EtaN (DYNAMICS)',myThid)
649           CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (DYNAMICS)',myThid)
650         CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (DYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (DYNAMICS)',myThid)
651         CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (DYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (DYNAMICS)',myThid)
652         CALL DEBUG_STATS_RL(Nr,theta,'Theta (DYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,theta,'Theta (DYNAMICS)',myThid)
653         CALL DEBUG_STATS_RL(Nr,salt,'Salt (DYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,salt,'Salt (DYNAMICS)',myThid)
654         CALL DEBUG_STATS_RL(Nr,Gu,'Gu (DYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,gU,'Gu (DYNAMICS)',myThid)
655         CALL DEBUG_STATS_RL(Nr,Gv,'Gv (DYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,gV,'Gv (DYNAMICS)',myThid)
656         CALL DEBUG_STATS_RL(Nr,Gt,'Gt (DYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,gT,'Gt (DYNAMICS)',myThid)
657         CALL DEBUG_STATS_RL(Nr,Gs,'Gs (DYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,gS,'Gs (DYNAMICS)',myThid)
658         CALL DEBUG_STATS_RL(Nr,GuNm1,'GuNm1 (DYNAMICS)',myThid)  #ifndef ALLOW_ADAMSBASHFORTH_3
659         CALL DEBUG_STATS_RL(Nr,GvNm1,'GvNm1 (DYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,guNm1,'GuNm1 (DYNAMICS)',myThid)
660         CALL DEBUG_STATS_RL(Nr,GtNm1,'GtNm1 (DYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,gvNm1,'GvNm1 (DYNAMICS)',myThid)
661         CALL DEBUG_STATS_RL(Nr,GsNm1,'GsNm1 (DYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,gtNm1,'GtNm1 (DYNAMICS)',myThid)
662           CALL DEBUG_STATS_RL(Nr,gsNm1,'GsNm1 (DYNAMICS)',myThid)
663    #endif
664          ENDIF
665    #endif
666    
667    #ifdef DYNAMICS_GUGV_EXCH_CHECK
668    C- jmc: For safety checking only: This Exchange here should not change
669    C       the solution. If solution changes, it means something is wrong,
670    C       but it does not mean that it is less wrong with this exchange.
671          IF ( debugLevel .GT. debLevB ) THEN
672           CALL EXCH_UV_XYZ_RL(gU,gV,.TRUE.,myThid)
673          ENDIF
674    #endif
675    
676    #ifdef ALLOW_DEBUG
677          IF ( debugLevel .GE. debLevB )
678         &   CALL DEBUG_LEAVE( 'DYNAMICS', myThid )
679  #endif  #endif
680    
681        RETURN        RETURN
Legend:
Removed from v.1.69  
changed lines
  Added in v.1.133
  ViewVC Help
Powered by ViewVC 1.1.22