/[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.167 by jmc, Tue Nov 5 13:34:31 2013 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_MOM_COMMON
7    # include "MOM_COMMON_OPTIONS.h"
8    #endif
9    #ifdef ALLOW_OBCS
10    # include "OBCS_OPTIONS.h"
11    #endif
12    
13    #undef DYNAMICS_GUGV_EXCH_CHECK
14    
15    CBOP
16    C     !ROUTINE: DYNAMICS
17    C     !INTERFACE:
18        SUBROUTINE DYNAMICS(myTime, myIter, myThid)        SUBROUTINE DYNAMICS(myTime, myIter, myThid)
19  C     /==========================================================\  C     !DESCRIPTION: \bv
20  C     | SUBROUTINE DYNAMICS                                      |  C     *==========================================================*
21  C     | o Controlling routine for the explicit part of the model |  C     | SUBROUTINE DYNAMICS
22  C     |   dynamics.                                              |  C     | o Controlling routine for the explicit part of the model
23  C     |==========================================================|  C     |   dynamics.
24  C     | This routine evaluates the "dynamics" terms for each     |  C     *==========================================================*
25  C     | block of ocean in turn. Because the blocks of ocean have |  C     | This routine evaluates the "dynamics" terms for each
26  C     | overlap regions they are independent of one another.     |  C     | block of ocean in turn. Because the blocks of ocean have
27  C     | If terms involving lateral integrals are needed in this  |  C     | overlap regions they are independent of one another.
28  C     | routine care will be needed. Similarly finite-difference |  C     | If terms involving lateral integrals are needed in this
29  C     | operations with stencils wider than the overlap region   |  C     | routine care will be needed. Similarly finite-difference
30  C     | require special consideration.                           |  C     | operations with stencils wider than the overlap region
31  C     | Notes                                                    |  C     | require special consideration.
32  C     | =====                                                    |  C     | The algorithm...
33  C     | C*P* comments indicating place holders for which code is |  C     |
34  C     |      presently being developed.                          |  C     | "Correction Step"
35  C     \==========================================================/  C     | =================
36    C     | Here we update the horizontal velocities with the surface
37    C     | pressure such that the resulting flow is either consistent
38    C     | with the free-surface evolution or the rigid-lid:
39    C     |   U[n] = U* + dt x d/dx P
40    C     |   V[n] = V* + dt x d/dy P
41    C     |   W[n] = W* + dt x d/dz P  (NH mode)
42    C     |
43    C     | "Calculation of Gs"
44    C     | ===================
45    C     | This is where all the accelerations and tendencies (ie.
46    C     | physics, parameterizations etc...) are calculated
47    C     |   rho = rho ( theta[n], salt[n] )
48    C     |   b   = b(rho, theta)
49    C     |   K31 = K31 ( rho )
50    C     |   Gu[n] = Gu( u[n], v[n], wVel, b, ... )
51    C     |   Gv[n] = Gv( u[n], v[n], wVel, b, ... )
52    C     |   Gt[n] = Gt( theta[n], u[n], v[n], wVel, K31, ... )
53    C     |   Gs[n] = Gs( salt[n], u[n], v[n], wVel, K31, ... )
54    C     |
55    C     | "Time-stepping" or "Prediction"
56    C     | ================================
57    C     | The models variables are stepped forward with the appropriate
58    C     | time-stepping scheme (currently we use Adams-Bashforth II)
59    C     | - For momentum, the result is always *only* a "prediction"
60    C     | in that the flow may be divergent and will be "corrected"
61    C     | later with a surface pressure gradient.
62    C     | - Normally for tracers the result is the new field at time
63    C     | level [n+1} *BUT* in the case of implicit diffusion the result
64    C     | is also *only* a prediction.
65    C     | - We denote "predictors" with an asterisk (*).
66    C     |   U* = U[n] + dt x ( 3/2 Gu[n] - 1/2 Gu[n-1] )
67    C     |   V* = V[n] + dt x ( 3/2 Gv[n] - 1/2 Gv[n-1] )
68    C     |   theta[n+1] = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
69    C     |   salt[n+1] = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
70    C     | With implicit diffusion:
71    C     |   theta* = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
72    C     |   salt* = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
73    C     |   (1 + dt * K * d_zz) theta[n] = theta*
74    C     |   (1 + dt * K * d_zz) salt[n] = salt*
75    C     |
76    C     *==========================================================*
77    C     \ev
78    C     !USES:
79        IMPLICIT NONE        IMPLICIT NONE
   
80  C     == Global variables ===  C     == Global variables ===
81  #include "SIZE.h"  #include "SIZE.h"
82  #include "EEPARAMS.h"  #include "EEPARAMS.h"
83  #include "PARAMS.h"  #include "PARAMS.h"
 #include "DYNVARS.h"  
84  #include "GRID.h"  #include "GRID.h"
85    #include "DYNVARS.h"
86    #ifdef ALLOW_MOM_COMMON
87    # include "MOM_VISC.h"
88    #endif
89    #ifdef ALLOW_CD_CODE
90    # include "CD_CODE_VARS.h"
91    #endif
92  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
93  # include "tamc.h"  # include "tamc.h"
94  # include "tamc_keys.h"  # include "tamc_keys.h"
95  # include "FFIELDS.h"  # include "FFIELDS.h"
96    # include "EOS.h"
97  # ifdef ALLOW_KPP  # ifdef ALLOW_KPP
98  #  include "KPP.h"  #  include "KPP.h"
99  # endif  # endif
100  # ifdef ALLOW_GMREDI  # ifdef ALLOW_PTRACERS
101  #  include "GMREDI.h"  #  include "PTRACERS_SIZE.h"
102    #  include "PTRACERS_FIELDS.h"
103    # endif
104    # ifdef ALLOW_OBCS
105    #  include "OBCS_PARAMS.h"
106    #  include "OBCS_FIELDS.h"
107    #  ifdef ALLOW_PTRACERS
108    #   include "OBCS_PTRACERS.h"
109    #  endif
110    # endif
111    # ifdef ALLOW_MOM_FLUXFORM
112    #  include "MOM_FLUXFORM.h"
113  # endif  # endif
114  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
115    
116  #ifdef ALLOW_TIMEAVE  C     !CALLING SEQUENCE:
117  #include "TIMEAVE_STATV.h"  C     DYNAMICS()
118  #endif  C      |
119    C      |-- CALC_EP_FORCING
120    C      |
121    C      |-- CALC_GRAD_PHI_SURF
122    C      |
123    C      |-- CALC_VISCOSITY
124    C      |
125    C      |-- MOM_CALC_3D_STRAIN
126    C      |
127    C      |-- CALC_EDDY_STRESS
128    C      |
129    C      |-- CALC_PHI_HYD
130    C      |
131    C      |-- MOM_FLUXFORM
132    C      |
133    C      |-- MOM_VECINV
134    C      |
135    C      |-- MOM_CALC_SMAG_3D
136    C      |-- MOM_UV_SMAG_3D
137    C      |
138    C      |-- TIMESTEP
139    C      |
140    C      |-- MOM_U_IMPLICIT_R
141    C      |-- MOM_V_IMPLICIT_R
142    C      |
143    C      |-- IMPLDIFF
144    C      |
145    C      |-- OBCS_APPLY_UV
146    C      |
147    C      |-- CALC_GW
148    C      |
149    C      |-- DIAGNOSTICS_FILL
150    C      |-- DEBUG_STATS_RL
151    
152    C     !INPUT/OUTPUT PARAMETERS:
153  C     == Routine arguments ==  C     == Routine arguments ==
154  C     myTime - Current time in simulation  C     myTime :: Current time in simulation
155  C     myIter - Current iteration number in simulation  C     myIter :: Current iteration number in simulation
156  C     myThid - Thread number for this instance of the routine.  C     myThid :: Thread number for this instance of the routine.
157        _RL myTime        _RL myTime
158        INTEGER myIter        INTEGER myIter
159        INTEGER myThid        INTEGER myThid
160    
161    C     !FUNCTIONS:
162    #ifdef ALLOW_DIAGNOSTICS
163          LOGICAL  DIAGNOSTICS_IS_ON
164          EXTERNAL DIAGNOSTICS_IS_ON
165    #endif
166    
167    C     !LOCAL VARIABLES:
168  C     == Local variables  C     == Local variables
169  C     xA, yA                 - Per block temporaries holding face areas  C     fVer[UV]               o fVer: Vertical flux term - note fVer
170  C     uTrans, vTrans, rTrans - Per block temporaries holding flow  C                                    is "pipelined" in the vertical
171  C                              transport  C                                    so we need an fVer for each
172  C                              o uTrans: Zonal transport  C                                    variable.
173  C                              o vTrans: Meridional transport  C     phiHydC    :: hydrostatic potential anomaly at cell center
174  C                              o rTrans: Vertical transport  C                   In z coords phiHyd is the hydrostatic potential
175  C     maskUp                   o maskUp: land/water mask for W points  C                      (=pressure/rho0) anomaly
176  C     fVer[STUV]               o fVer: Vertical flux term - note fVer  C                   In p coords phiHyd is the geopotential height anomaly.
177  C                                      is "pipelined" in the vertical  C     phiHydF    :: hydrostatic potential anomaly at middle between 2 centers
178  C                                      so we need an fVer for each  C     dPhiHydX,Y :: Gradient (X & Y directions) of hydrostatic potential anom.
179  C                                      variable.  C     phiSurfX,  ::  gradient of Surface potential (Pressure/rho, ocean)
180  C     rhoK, rhoKM1   - Density at current level, and level above  C     phiSurfY             or geopotential (atmos) in X and Y direction
181  C     phiHyd         - Hydrostatic part of the potential phiHydi.  C     guDissip   :: dissipation tendency (all explicit terms), u component
182  C                      In z coords phiHydiHyd is the hydrostatic  C     gvDissip   :: dissipation tendency (all explicit terms), v component
183  C                      Potential (=pressure/rho0) anomaly  C     KappaRU    :: vertical viscosity for velocity U-component
184  C                      In p coords phiHydiHyd is the geopotential  C     KappaRV    :: vertical viscosity for velocity V-component
185  C                      surface height anomaly.  C     iMin, iMax :: Ranges and sub-block indices on which calculations
186  C     phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean)  C     jMin, jMax    are applied.
187  C     phiSurfY             or geopotentiel (atmos) in X and Y direction  C     bi, bj     :: tile indices
188  C     KappaRT,       - Total diffusion in vertical for T and S.  C     k          :: current level index
189  C     KappaRS          (background + spatially varying, isopycnal term).  C     km1, kp1   :: index of level above (k-1) and below (k+1)
190  C     iMin, iMax     - Ranges and sub-block indices on which calculations  C     kUp, kDown :: Index for interface above and below. kUp and kDown are
191  C     jMin, jMax       are applied.  C                   are switched with k to be the appropriate index into fVerU,V
 C     bi, bj  
 C     k, kup,        - Index for layer above and below. kup and kDown  
 C     kDown, km1       are switched with layer to be the appropriate  
 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)  
192        _RL fVerU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
193        _RL fVerV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
194        _RL phiHyd  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL phiHydF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
195        _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiHydC (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
196        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL dPhiHydX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
197          _RL dPhiHydY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
198        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
199        _RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
200        _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL guDissip(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
201        _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL gvDissip(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
202        _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
203        _RL KappaRV (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
204        _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  #ifdef ALLOW_SMAG_3D
205        _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  C     str11       :: strain component Vxx @ grid-cell center
206        _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  C     str22       :: strain component Vyy @ grid-cell center
207        _RL tauAB  C     str33       :: strain component Vzz @ grid-cell center
208    C     str12       :: strain component Vxy @ grid-cell corner
209  C This is currently used by IVDC and Diagnostics  C     str13       :: strain component Vxz @ above uVel
210        _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  C     str23       :: strain component Vyz @ above vVel
211    C     viscAh3d_00 :: Smagorinsky viscosity @ grid-cell center
212    C     viscAh3d_12 :: Smagorinsky viscosity @ grid-cell corner
213    C     viscAh3d_13 :: Smagorinsky viscosity @ above uVel
214    C     viscAh3d_23 :: Smagorinsky viscosity @ above vVel
215    C     addDissU    :: zonal momentum tendency from 3-D Smag. viscosity
216    C     addDissV    :: merid momentum tendency from 3-D Smag. viscosity
217          _RL str11(1-OLx:sNx+OLx,1-OLy:sNy+OLy, Nr )
218          _RL str22(1-OLx:sNx+OLx,1-OLy:sNy+OLy, Nr )
219          _RL str33(1-OLx:sNx+OLx,1-OLy:sNy+OLy, Nr )
220          _RL str12(1-OLx:sNx+OLx,1-OLy:sNy+OLy, Nr )
221          _RL str13(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)
222          _RL str23(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)
223          _RL viscAh3d_00(1-OLx:sNx+OLx,1-OLy:sNy+OLy, Nr )
224          _RL viscAh3d_12(1-OLx:sNx+OLx,1-OLy:sNy+OLy, Nr )
225          _RL viscAh3d_13(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)
226          _RL viscAh3d_23(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)
227          _RL addDissU(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
228          _RL addDissV(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
229    #elif ( defined ALLOW_NONHYDROSTATIC )
230          _RL str13(1), str23(1), str33(1)
231          _RL viscAh3d_00(1), viscAh3d_13(1), viscAh3d_23(1)
232    #endif
233    
       INTEGER iMin, iMax  
       INTEGER jMin, jMax  
234        INTEGER bi, bj        INTEGER bi, bj
235        INTEGER i, j        INTEGER i, j
236        INTEGER k, km1, kup, kDown        INTEGER k, km1, kp1, kUp, kDown
237          INTEGER iMin, iMax
238          INTEGER jMin, jMax
239          PARAMETER( iMin = 0 , iMax = sNx+1 )
240          PARAMETER( jMin = 0 , jMax = sNy+1 )
241    
242    #ifdef ALLOW_DIAGNOSTICS
243          LOGICAL dPhiHydDiagIsOn
244          _RL tmpFac
245    #endif /* ALLOW_DIAGNOSTICS */
246    
 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)  
   
247  C---    The algorithm...  C---    The algorithm...
248  C  C
249  C       "Correction Step"  C       "Correction Step"
# Line 165  C         salt* = salt[n] + dt x ( 3/2 G Line 287  C         salt* = salt[n] + dt x ( 3/2 G
287  C         (1 + dt * K * d_zz) theta[n] = theta*  C         (1 + dt * K * d_zz) theta[n] = theta*
288  C         (1 + dt * K * d_zz) salt[n] = salt*  C         (1 + dt * K * d_zz) salt[n] = salt*
289  C---  C---
290    CEOP
291    
292  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_DEBUG
293  C--   dummy statement to end declaration part        IF (debugMode) CALL DEBUG_ENTER( 'DYNAMICS', myThid )
294        ikey = 1  #endif
 #endif /* ALLOW_AUTODIFF_TAMC */  
295    
296  C--   Set up work arrays with valid (i.e. not NaN) values  #ifdef ALLOW_DIAGNOSTICS
297  C     These inital values do not alter the numerical results. They        dPhiHydDiagIsOn = .FALSE.
298  C     just ensure that all memory references are to valid floating        IF ( useDiagnostics )
299  C     point numbers. This prevents spurious hardware signals due to       &  dPhiHydDiagIsOn = DIAGNOSTICS_IS_ON( 'Um_dPHdx', myThid )
300  C     uninitialised but inert locations.       &               .OR. DIAGNOSTICS_IS_ON( 'Vm_dPHdy', myThid )
301        DO j=1-OLy,sNy+OLy  #endif
302         DO i=1-OLx,sNx+OLx  
303          xA(i,j)      = 0. _d 0  C-- Call to routine for calculation of Eliassen-Palm-flux-forced
304          yA(i,j)      = 0. _d 0  C    U-tendency, if desired:
305          uTrans(i,j)  = 0. _d 0  #ifdef INCLUDE_EP_FORCING_CODE
306          vTrans(i,j)  = 0. _d 0        CALL CALC_EP_FORCING(myThid)
307          DO k=1,Nr  #endif
          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  
308    
309    #ifdef ALLOW_AUTODIFF_MONITOR_DIAG
310          CALL DUMMY_IN_DYNAMICS( myTime, myIter, myThid )
311    #endif
312    
313  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
314  C--   HPF directive to help TAMC  C--   HPF directive to help TAMC
# Line 207  CHPF$ INDEPENDENT Line 319  CHPF$ INDEPENDENT
319    
320  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
321  C--    HPF directive to help TAMC  C--    HPF directive to help TAMC
322  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS,fVerU,fVerV  CHPF$  INDEPENDENT, NEW (fVerU,fVerV
323  CHPF$&                  ,phiHyd,utrans,vtrans,xA,yA  CHPF$&                  ,phiHydF
324  CHPF$&                  ,KappaRT,KappaRS,KappaRU,KappaRV  CHPF$&                  ,KappaRU,KappaRV
325  CHPF$&                  )  CHPF$&                  )
326  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
327    
# Line 218  CHPF$&                  ) Line 330  CHPF$&                  )
330  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
331            act1 = bi - myBxLo(myThid)            act1 = bi - myBxLo(myThid)
332            max1 = myBxHi(myThid) - myBxLo(myThid) + 1            max1 = myBxHi(myThid) - myBxLo(myThid) + 1
   
333            act2 = bj - myByLo(myThid)            act2 = bj - myByLo(myThid)
334            max2 = myByHi(myThid) - myByLo(myThid) + 1            max2 = myByHi(myThid) - myByLo(myThid) + 1
   
335            act3 = myThid - 1            act3 = myThid - 1
336            max3 = nTx*nTy            max3 = nTx*nTy
   
337            act4 = ikey_dynamics - 1            act4 = ikey_dynamics - 1
338              idynkey = (act1 + 1) + act2*max1
           ikey = (act1 + 1) + act2*max1  
339       &                      + act3*max1*max2       &                      + act3*max1*max2
340       &                      + act4*max1*max2*max3       &                      + act4*max1*max2*max3
341  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
342    
343  C--     Set up work arrays that need valid initial values  C--   Set up work arrays with valid (i.e. not NaN) values
344          DO j=1-OLy,sNy+OLy  C     These initial values do not alter the numerical results. They
345           DO i=1-OLx,sNx+OLx  C     just ensure that all memory references are to valid floating
346            rTrans(i,j)   = 0. _d 0  C     point numbers. This prevents spurious hardware signals due to
347            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  
348    
349    #ifdef ALLOW_AUTODIFF_TAMC
350          DO k=1,Nr          DO k=1,Nr
351           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
352            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
353  C This is currently also used by IVDC and Diagnostics  cph(
354             ConvectCount(i,j,k) = 0.  c--   need some re-initialisation here to break dependencies
355             KappaRT(i,j,k) = 0. _d 0  cph)
356             KappaRS(i,j,k) = 0. _d 0             gU(i,j,k,bi,bj) = 0. _d 0
357               gV(i,j,k,bi,bj) = 0. _d 0
358            ENDDO            ENDDO
359           ENDDO           ENDDO
360          ENDDO          ENDDO
   
         iMin = 1-OLx+1  
         iMax = sNx+OLx  
         jMin = 1-OLy+1  
         jMax = sNy+OLy  
   
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
 CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
 CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
 CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 C--     Start of diagnostic loop  
         DO k=Nr,1,-1  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 C? Patrick, is this formula correct now that we change the loop range?  
 C? Do we still need this?  
 cph kkey formula corrected.  
 cph Needed for rhok, rhokm1, in the case useGMREDI.  
          kkey = (ikey-1)*Nr + k  
 CADJ STORE rhokm1(:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte  
 CADJ STORE rhok  (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 C--       Integrate continuity vertically for vertical velocity  
           CALL INTEGRATE_FOR_W(  
      I                         bi, bj, k, uVel, vVel,  
      O                         wVel,  
      I                         myThid )  
   
 #ifdef    ALLOW_OBCS  
 #ifdef    ALLOW_NONHYDROSTATIC  
 C--       Apply OBC to W if in N-H mode  
           IF (useOBCS.AND.nonHydrostatic) THEN  
             CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )  
           ENDIF  
 #endif    /* ALLOW_NONHYDROSTATIC */  
 #endif    /* ALLOW_OBCS */  
   
 C--       Calculate gradients of potential density for isoneutral  
 C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)  
 c         IF ( k.GT.1 .AND. (useGMRedi.OR.ivdc_kappa.NE.0.) ) THEN  
           IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.) ) THEN  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  
 CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
             CALL FIND_RHO(  
      I        bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,  
      I        theta, salt,  
      O        rhoK,  
      I        myThid )  
             IF (k.GT.1) THEN  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  
 CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  
361  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
362               CALL FIND_RHO(          DO j=1-OLy,sNy+OLy
363       I        bi, bj, iMin, iMax, jMin, jMax, k-1, k, eosType,           DO i=1-OLx,sNx+OLx
364       I        theta, salt,            fVerU  (i,j,1) = 0. _d 0
365       O        rhoKm1,            fVerU  (i,j,2) = 0. _d 0
366       I        myThid )            fVerV  (i,j,1) = 0. _d 0
367              ENDIF            fVerV  (i,j,2) = 0. _d 0
368              CALL GRAD_SIGMA(            phiHydF (i,j)  = 0. _d 0
369       I             bi, bj, iMin, iMax, jMin, jMax, k,            phiHydC (i,j)  = 0. _d 0
370       I             rhoK, rhoKm1, rhoK,  #ifndef INCLUDE_PHIHYD_CALCULATION_CODE
371       O             sigmaX, sigmaY, sigmaR,            dPhiHydX(i,j)  = 0. _d 0
372       I             myThid )            dPhiHydY(i,j)  = 0. _d 0
373            ENDIF  #endif
374              phiSurfX(i,j)  = 0. _d 0
375  C--       Implicit Vertical Diffusion for Convection            phiSurfY(i,j)  = 0. _d 0
376  c ==> should use sigmaR !!!            guDissip(i,j)  = 0. _d 0
377            IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN            gvDissip(i,j)  = 0. _d 0
378              CALL CALC_IVDC(  #ifdef ALLOW_AUTODIFF_TAMC
379       I        bi, bj, iMin, iMax, jMin, jMax, k,            phiHydLow(i,j,bi,bj) = 0. _d 0
380       I        rhoKm1, rhoK,  # if (defined NONLIN_FRSURF) && (defined ALLOW_MOM_FLUXFORM)
381       U        ConvectCount, KappaRT, KappaRS,  #  ifndef DISABLE_RSTAR_CODE
382       I        myTime, myIter, myThid)  #   ifndef ALLOW_AUTODIFF_OPENAD
383            ENDIF            dWtransC(i,j,bi,bj) = 0. _d 0
384              dWtransU(i,j,bi,bj) = 0. _d 0
385  C--     end of diagnostic k loop (Nr:1)            dWtransV(i,j,bi,bj) = 0. _d 0
386    #   endif
387    #  endif
388    # endif
389    #endif
390             ENDDO
391          ENDDO          ENDDO
392    
393    C--     Start computation of dynamics
394    
395  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
396  cph avoids recomputation of integrate_for_w  CADJ STORE wVel (:,:,:,bi,bj) =
397  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ &     comlev1_bibj, key=idynkey, byte=isbyte
398  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
399    
400  #ifdef  ALLOW_OBCS  C--     Explicit part of the Surface Potential Gradient (add in TIMESTEP)
401  C--     Calculate future values on open boundaries  C       (note: this loop will be replaced by CALL CALC_GRAD_ETA)
402          IF (useOBCS) THEN          IF (implicSurfPress.NE.1.) THEN
403            CALL OBCS_CALC( bi, bj, myTime+deltaT,            CALL CALC_GRAD_PHI_SURF(
404       I            uVel, vVel, wVel, theta, salt,       I         bi,bj,iMin,iMax,jMin,jMax,
405       I            myThid )       I         etaN,
406         O         phiSurfX,phiSurfY,
407         I         myThid )
408          ENDIF          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  
409    
410  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
411  CADJ STORE sigmaX(:,:,:) = comlev1, key=ikey, byte=isbyte  CADJ STORE uVel (:,:,:,bi,bj) = comlev1_bibj, key=idynkey, byte=isbyte
412  CADJ STORE sigmaY(:,:,:) = comlev1, key=ikey, byte=isbyte  CADJ STORE vVel (:,:,:,bi,bj) = comlev1_bibj, key=idynkey, byte=isbyte
413  CADJ STORE sigmaR(:,:,:) = comlev1, key=ikey, byte=isbyte  #ifdef ALLOW_KPP
414    CADJ STORE KPPviscAz (:,:,:,bi,bj)
415    CADJ &                 = comlev1_bibj, key=idynkey, byte=isbyte
416    #endif /* ALLOW_KPP */
417  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
418  C--     Calculate iso-neutral slopes for the GM/Redi parameterisation  
419          IF (useGMRedi) THEN  #if (defined INCLUDE_CALC_DIFFUSIVITY_CALL) && !(defined ALLOW_AUTODIFF)
420            IF ( .NOT.momViscosity ) THEN
421    #endif /* INCLUDE_CALC_DIFFUSIVITY_CALL and not ALLOW_AUTODIFF */
422            DO k=1,Nr            DO k=1,Nr
423              CALL GMREDI_CALC_TENSOR(             DO j=1-OLy,sNy+OLy
424       I             bi, bj, iMin, iMax, jMin, jMax, k,              DO i=1-OLx,sNx+OLx
425       I             sigmaX, sigmaY, sigmaR,               KappaRU(i,j,k) = 0. _d 0
426       I             myThid )               KappaRV(i,j,k) = 0. _d 0
427                ENDDO
428               ENDDO
429            ENDDO            ENDDO
430  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef INCLUDE_CALC_DIFFUSIVITY_CALL
431    C--     Calculate the total vertical viscosity
432    #ifdef ALLOW_AUTODIFF
433            IF ( momViscosity ) THEN
434    #else
435          ELSE          ELSE
436            DO k=1, Nr  #endif
437              CALL GMREDI_CALC_TENSOR_DUMMY(            CALL CALC_VISCOSITY(
438       I             bi, bj, iMin, iMax, jMin, jMax, k,       I            bi,bj, iMin,iMax,jMin,jMax,
439       I             sigmaX, sigmaY, sigmaR,       O            KappaRU, KappaRV,
440       I             myThid )       I            myThid )
           ENDDO  
 #endif /* ALLOW_AUTODIFF_TAMC */  
441          ENDIF          ENDIF
442    #endif /* INCLUDE_CALC_DIFFUSIVITY_CALL */
443    
444  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_SMAG_3D
445  CADJ STORE Kwx(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte          IF ( useSmag3D ) THEN
446  CADJ STORE Kwy(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte            CALL MOM_CALC_3D_STRAIN(
447  CADJ STORE Kwz(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte       O         str11, str22, str33, str12, str13, str23,
448  #endif /* ALLOW_AUTODIFF_TAMC */       I         bi, bj, myThid )
   
 #endif  /* ALLOW_GMREDI */  
   
 #ifdef  ALLOW_KPP  
 C--     Compute KPP mixing coefficients  
         IF (useKPP) THEN  
           CALL KPP_CALC(  
      I                  bi, bj, myTime, myThid )  
 #ifdef ALLOW_AUTODIFF_TAMC  
         ELSE  
           CALL KPP_CALC_DUMMY(  
      I                  bi, bj, myTime, myThid )  
 #endif /* ALLOW_AUTODIFF_TAMC */  
449          ENDIF          ENDIF
450    #endif /* ALLOW_SMAG_3D */
451    
452  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
453  CADJ STORE KPPghat   (:,:,:,bi,bj)  CADJ STORE KappaRU(:,:,:)
454  CADJ &   , KPPviscAz (:,:,:,bi,bj)  CADJ &     = comlev1_bibj, key=idynkey, byte=isbyte
455  CADJ &   , KPPdiffKzT(:,:,:,bi,bj)  CADJ STORE KappaRV(:,:,:)
456  CADJ &   , KPPdiffKzS(:,:,:,bi,bj)  CADJ &     = comlev1_bibj, key=idynkey, byte=isbyte
457  CADJ &   , KPPfrac   (:,:  ,bi,bj)  #endif /* ALLOW_AUTODIFF_TAMC */
458  CADJ &                 = comlev1_bibj, key=ikey, byte=isbyte  
459  #endif /* ALLOW_AUTODIFF_TAMC */  #ifdef ALLOW_OBCS
460    C--   For Stevens boundary conditions velocities need to be extrapolated
461  #endif  /* ALLOW_KPP */  C     (copied) to a narrow strip outside the domain
462            IF ( useOBCS ) THEN
463  #ifdef ALLOW_AUTODIFF_TAMC            CALL OBCS_COPY_UV_N(
464  CADJ STORE KappaRT(:,:,:)     = comlev1_bibj, key = ikey, byte = isbyte       U         uVel(1-OLx,1-OLy,1,bi,bj),
465  CADJ STORE KappaRS(:,:,:)     = comlev1_bibj, key = ikey, byte = isbyte       U         vVel(1-OLx,1-OLy,1,bi,bj),
466  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte       I         Nr, bi, bj, myThid )
 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)  
467          ENDIF          ENDIF
468  #endif /* ALLOW_AIM */  #endif /* ALLOW_OBCS */
469    
470    #ifdef ALLOW_EDDYPSI
471            CALL CALC_EDDY_STRESS(bi,bj,myThid)
472    #endif
473    
474  C--     Start of thermodynamics loop  C--     Start of dynamics loop
475          DO k=Nr,1,-1          DO k=1,Nr
 #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 */  
476    
477  C--       km1    Points to level above k (=k-1)  C--       km1    Points to level above k (=k-1)
478  C--       kup    Cycles through 1,2 to point to layer above  C--       kup    Cycles through 1,2 to point to layer above
479  C--       kDown  Cycles through 2,1 to point to current layer  C--       kDown  Cycles through 2,1 to point to current layer
480    
481            km1  = MAX(1,k-1)            km1  = MAX(1,k-1)
482              kp1  = MIN(k+1,Nr)
483            kup  = 1+MOD(k+1,2)            kup  = 1+MOD(k+1,2)
484            kDown= 1+MOD(k,2)            kDown= 1+MOD(k,2)
485    
           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  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL  
 C--      Calculate the total vertical diffusivity  
          CALL CALC_DIFFUSIVITY(  
      I        bi,bj,iMin,iMax,jMin,jMax,k,  
      I        maskUp,  
      O        KappaRT,KappaRS,KappaRU,KappaRV,  
      I        myThid)  
 #endif  
   
 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  
486  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
487  CADJ STORE gTNm1(:,:,k,bi,bj) = comlev1_bibj_k           kkey = (idynkey-1)*Nr + k
488  CADJ &   , key = kkey, byte = isbyte  c
489  #endif /* ALLOW_AUTODIFF_TAMC */  CADJ STORE totPhiHyd (:,:,k,bi,bj)
490              CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )  CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
491           END IF  CADJ STORE phiHydLow (:,:,bi,bj)
492    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
493  C--     end of thermodynamic k loop (Nr:1)  CADJ STORE theta (:,:,k,bi,bj)
494          ENDDO  CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
495    CADJ STORE salt  (:,:,k,bi,bj)
496    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
497  #ifdef ALLOW_AUTODIFF_TAMC  CADJ STORE gT(:,:,k,bi,bj)
498  C? Patrick? What about this one?  CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
499  cph Keys iikey and idkey don't seem to be needed  CADJ STORE gS(:,:,k,bi,bj)
500  cph since storing occurs on different tape for each  CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
501  cph impldiff call anyways.  # ifdef NONLIN_FRSURF
502  cph Thus, common block comlev1_impl isn't needed either.  cph-test
503  cph Storing below needed in the case useGMREDI.  CADJ STORE  phiHydC (:,:)
504          iikey = (ikey-1)*maximpl  CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
505  #endif /* ALLOW_AUTODIFF_TAMC */  CADJ STORE  phiHydF (:,:)
506    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
507  C--     Implicit diffusion  CADJ STORE  guDissip (:,:)
508          IF (implicitDiffusion) THEN  CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
509    CADJ STORE  gvDissip (:,:)
510           IF (tempStepping) THEN  CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
511  #ifdef ALLOW_AUTODIFF_TAMC  CADJ STORE  fVerU (:,:,:)
512              idkey = iikey + 1  CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
513  CADJ STORE gTNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE  fVerV (:,:,:)
514  #endif /* ALLOW_AUTODIFF_TAMC */  CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
515              CALL IMPLDIFF(  CADJ STORE gU(:,:,k,bi,bj)
516       I         bi, bj, iMin, iMax, jMin, jMax,  CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
517       I         deltaTtracer, KappaRT, recip_HFacC,  CADJ STORE gV(:,:,k,bi,bj)
518       U         gTNm1,  CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
519       I         myThid )  #  ifndef ALLOW_ADAMSBASHFORTH_3
520           ENDIF  CADJ STORE guNm1(:,:,k,bi,bj)
521    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
522           IF (saltStepping) THEN  CADJ STORE gvNm1(:,:,k,bi,bj)
523  #ifdef ALLOW_AUTODIFF_TAMC  CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
524           idkey = iikey + 2  #  else
525  CADJ STORE gSNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE guNm(:,:,k,bi,bj,1)
526    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
527    CADJ STORE guNm(:,:,k,bi,bj,2)
528    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
529    CADJ STORE gvNm(:,:,k,bi,bj,1)
530    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
531    CADJ STORE gvNm(:,:,k,bi,bj,2)
532    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
533    #  endif
534    #  ifdef ALLOW_CD_CODE
535    CADJ STORE uNM1(:,:,k,bi,bj)
536    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
537    CADJ STORE vNM1(:,:,k,bi,bj)
538    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
539    CADJ STORE uVelD(:,:,k,bi,bj)
540    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
541    CADJ STORE vVelD(:,:,k,bi,bj)
542    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
543    #  endif
544    # endif
545    # ifdef ALLOW_DEPTH_CONTROL
546    CADJ STORE  fVerU (:,:,:)
547    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
548    CADJ STORE  fVerV (:,:,:)
549    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
550    # endif
551  #endif /* ALLOW_AUTODIFF_TAMC */  #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  
552    
553  C--     Explicit part of the Surface Potentiel Gradient (add in TIMESTEP)  C--      Integrate hydrostatic balance for phiHyd with BC of phiHyd(z=0)=0
554  C       (note: this loop will be replaced by CALL CALC_GRAD_ETA)           IF ( implicitIntGravWave ) THEN
         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  
   
 C--     Start of dynamics loop  
         DO k=1,Nr  
   
 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)  
   
 C--      Integrate hydrostatic balance for phiHyd with BC of  
 C        phiHyd(z=0)=0  
 C        distinguishe between Stagger and Non Stagger time stepping  
          IF (staggerTimeStep) THEN  
555             CALL CALC_PHI_HYD(             CALL CALC_PHI_HYD(
556       I        bi,bj,iMin,iMax,jMin,jMax,k,       I        bi,bj,iMin,iMax,jMin,jMax,k,
557       I        gTnm1, gSnm1,       I        gT, gS,
558       U        phiHyd,       U        phiHydF,
559       I        myThid )       O        phiHydC, dPhiHydX, dPhiHydY,
560         I        myTime, myIter, myThid )
561           ELSE           ELSE
562             CALL CALC_PHI_HYD(             CALL CALC_PHI_HYD(
563       I        bi,bj,iMin,iMax,jMin,jMax,k,       I        bi,bj,iMin,iMax,jMin,jMax,k,
564       I        theta, salt,       I        theta, salt,
565       U        phiHyd,       U        phiHydF,
566       I        myThid )       O        phiHydC, dPhiHydX, dPhiHydY,
567         I        myTime, myIter, myThid )
568           ENDIF           ENDIF
569    #ifdef ALLOW_DIAGNOSTICS
570             IF ( dPhiHydDiagIsOn ) THEN
571               tmpFac = -1. _d 0
572               CALL DIAGNOSTICS_SCALE_FILL( dPhiHydX, tmpFac, 1,
573         &                           'Um_dPHdx', k, 1, 2, bi, bj, myThid )
574               CALL DIAGNOSTICS_SCALE_FILL( dPhiHydY, tmpFac, 1,
575         &                           'Vm_dPHdy', k, 1, 2, bi, bj, myThid )
576             ENDIF
577    #endif /* ALLOW_DIAGNOSTICS */
578    
579  C--      Calculate accelerations in the momentum equations (gU, gV, ...)  C--      Calculate accelerations in the momentum equations (gU, gV, ...)
580  C        and step forward storing the result in gUnm1, gVnm1, etc...  C        and step forward storing the result in gU, gV, etc...
581           IF ( momStepping ) THEN           IF ( momStepping ) THEN
582             CALL CALC_MOM_RHS(  #ifdef ALLOW_AUTODIFF_TAMC
583       I         bi,bj,iMin,iMax,jMin,jMax,k,kup,kDown,  # ifdef NONLIN_FRSURF
584       I         phiHyd,KappaRU,KappaRV,  #  if (defined ALLOW_MOM_FLUXFORM) && !(defined DISABLE_RSTAR_CODE)
585       U         fVerU, fVerV,  CADJ STORE dWtransC(:,:,bi,bj)
586       I         myTime, myThid)  CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
587    CADJ STORE dWtransU(:,:,bi,bj)
588    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
589    CADJ STORE dWtransV(:,:,bi,bj)
590    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
591    #  endif
592    CADJ STORE fVerU(:,:,:)
593    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
594    CADJ STORE fVerV(:,:,:)
595    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
596    # endif /* NONLIN_FRSURF */
597    #endif /* ALLOW_AUTODIFF_TAMC */
598               IF (.NOT. vectorInvariantMomentum) THEN
599    #ifdef ALLOW_MOM_FLUXFORM
600                  CALL MOM_FLUXFORM(
601         I         bi,bj,k,iMin,iMax,jMin,jMax,
602         I         KappaRU, KappaRV,
603         U         fVerU(1-OLx,1-OLy,kUp),   fVerV(1-OLx,1-OLy,kUp),
604         O         fVerU(1-OLx,1-OLy,kDown), fVerV(1-OLx,1-OLy,kDown),
605         O         guDissip, gvDissip,
606         I         myTime, myIter, myThid)
607    #endif
608               ELSE
609    #ifdef ALLOW_MOM_VECINV
610                 CALL MOM_VECINV(
611         I         bi,bj,k,iMin,iMax,jMin,jMax,
612         I         KappaRU, KappaRV,
613         I         fVerU(1-OLx,1-OLy,kUp),   fVerV(1-OLx,1-OLy,kUp),
614         O         fVerU(1-OLx,1-OLy,kDown), fVerV(1-OLx,1-OLy,kDown),
615         O         guDissip, gvDissip,
616         I         myTime, myIter, myThid)
617    #endif
618               ENDIF
619    
620    #ifdef ALLOW_SMAG_3D
621               IF ( useSmag3D ) THEN
622                 CALL MOM_CALC_SMAG_3D(
623         I         str11, str22, str33, str12, str13, str23,
624         O         viscAh3d_00, viscAh3d_12, viscAh3d_13, viscAh3d_23,
625         I         smag3D_hLsC, smag3D_hLsW, smag3D_hLsS, smag3D_hLsZ,
626         I         k, bi, bj, myThid )
627                 CALL MOM_UV_SMAG_3D(
628         I         str11, str22, str12, str13, str23,
629         I         viscAh3d_00, viscAh3d_12, viscAh3d_13, viscAh3d_23,
630         O         addDissU, addDissV,
631         I         iMin,iMax,jMin,jMax, k, bi, bj, myThid )
632                 DO j= jMin,jMax
633                  DO i= iMin,iMax
634                   guDissip(i,j) = guDissip(i,j) + addDissU(i,j)
635                   gvDissip(i,j) = gvDissip(i,j) + addDissV(i,j)
636                  ENDDO
637                 ENDDO
638               ENDIF
639    #endif /* ALLOW_SMAG_3D */
640    
641             CALL TIMESTEP(             CALL TIMESTEP(
642       I         bi,bj,iMin,iMax,jMin,jMax,k,       I         bi,bj,iMin,iMax,jMin,jMax,k,
643       I         phiHyd, phiSurfX, phiSurfY,       I         dPhiHydX,dPhiHydY, phiSurfX, phiSurfY,
644       I         myIter, myThid)       I         guDissip, gvDissip,
645         I         myTime, myIter, myThid)
646    
 #ifdef   ALLOW_OBCS  
 C--      Apply open boundary conditions  
          IF (useOBCS) THEN  
            CALL OBCS_APPLY_UV( bi, bj, k, gUnm1, gVnm1, myThid )  
          END IF  
 #endif   /* ALLOW_OBCS */  
   
 #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 */  
647           ENDIF           ENDIF
648    
   
649  C--     end of dynamics k loop (1:Nr)  C--     end of dynamics k loop (1:Nr)
650          ENDDO          ENDDO
651    
652    C--     Implicit Vertical advection & viscosity
653    #if (defined (INCLUDE_IMPLVERTADV_CODE) && \
654  C--     Implicit viscosity       defined (ALLOW_MOM_COMMON) && !(defined ALLOW_AUTODIFF_TAMC))
655          IF (implicitViscosity.AND.momStepping) THEN          IF ( momImplVertAdv ) THEN
656              CALL MOM_U_IMPLICIT_R( kappaRU,
657         I                           bi, bj, myTime, myIter, myThid )
658              CALL MOM_V_IMPLICIT_R( kappaRV,
659         I                           bi, bj, myTime, myIter, myThid )
660            ELSEIF ( implicitViscosity ) THEN
661    #else /* INCLUDE_IMPLVERTADV_CODE */
662            IF     ( implicitViscosity ) THEN
663    #endif /* INCLUDE_IMPLVERTADV_CODE */
664  #ifdef    ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
665            idkey = iikey + 3  CADJ STORE KappaRU(:,:,:) = comlev1_bibj , key=idynkey, byte=isbyte
666  CADJ STORE gUNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE gU(:,:,:,bi,bj) = comlev1_bibj , key=idynkey, byte=isbyte
667  #endif    /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
668            CALL IMPLDIFF(            CALL IMPLDIFF(
669       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
670       I         deltaTmom, KappaRU,recip_HFacW,       I         -1, KappaRU, recip_hFacW(1-OLx,1-OLy,1,bi,bj),
671       U         gUNm1,       U         gU,
672       I         myThid )       I         myThid )
673  #ifdef    ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
674            idkey = iikey + 4  CADJ STORE KappaRV(:,:,:) = comlev1_bibj , key=idynkey, byte=isbyte
675  CADJ STORE gVNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE gV(:,:,:,bi,bj) = comlev1_bibj , key=idynkey, byte=isbyte
676  #endif    /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
677            CALL IMPLDIFF(            CALL IMPLDIFF(
678       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
679       I         deltaTmom, KappaRV,recip_HFacS,       I         -2, KappaRV, recip_hFacS(1-OLx,1-OLy,1,bi,bj),
680       U         gVNm1,       U         gV,
681       I         myThid )       I         myThid )
682            ENDIF
683    
684  #ifdef   ALLOW_OBCS  #ifdef ALLOW_OBCS
685  C--      Apply open boundary conditions  C--      Apply open boundary conditions
686           IF (useOBCS) THEN          IF ( useOBCS ) THEN
687             DO K=1,Nr  C--      but first save intermediate velocities to be used in the
688               CALL OBCS_APPLY_UV( bi, bj, k, gUnm1, gVnm1, myThid )  C        next time step for the Stevens boundary conditions
689             ENDDO            CALL OBCS_SAVE_UV_N(
690           END IF       I        bi, bj, iMin, iMax, jMin, jMax, 0,
691  #endif   /* ALLOW_OBCS */       I        gU, gV, myThid )
692              CALL OBCS_APPLY_UV( bi, bj, 0, gU, gV, myThid )
693            ENDIF
694    #endif /* ALLOW_OBCS */
695    
696  #ifdef    INCLUDE_CD_CODE  #ifdef    ALLOW_CD_CODE
697            IF (implicitViscosity.AND.useCDscheme) THEN
698  #ifdef    ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
699            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  
700  #endif    /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
701            CALL IMPLDIFF(            CALL IMPLDIFF(
702       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
703       I         deltaTmom, KappaRU,recip_HFacW,       I         0, KappaRU, recip_hFacW(1-OLx,1-OLy,1,bi,bj),
704       U         vVelD,       U         vVelD,
705       I         myThid )       I         myThid )
706  #ifdef    ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
707            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  
708  #endif    /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
709            CALL IMPLDIFF(            CALL IMPLDIFF(
710       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
711       I         deltaTmom, KappaRV,recip_HFacS,       I         0, KappaRV, recip_hFacS(1-OLx,1-OLy,1,bi,bj),
712       U         uVelD,       U         uVelD,
713       I         myThid )       I         myThid )
 #endif    /* INCLUDE_CD_CODE */  
 C--     End If implicitViscosity.AND.momStepping  
714          ENDIF          ENDIF
715    #endif    /* ALLOW_CD_CODE */
716  Cjmc : add for phiHyd output <- but not working if multi tile per CPU  C--     End implicit Vertical advection & viscosity
717  c       IF ( DIFFERENT_MULTIPLE(dumpFreq,myTime+deltaTClock,myTime)  
718  c    &  .AND. buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
719  c         WRITE(suff,'(I10.10)') myIter+1  
720  c         CALL WRITE_FLD_XYZ_RL('PH.',suff,phiHyd,myIter+1,myThid)  #ifdef ALLOW_NONHYDROSTATIC
721  c       ENDIF  C--   Step forward W field in N-H algorithm
722  Cjmc(end)          IF ( nonHydrostatic ) THEN
723    #ifdef ALLOW_DEBUG
724  #ifdef ALLOW_TIMEAVE           IF (debugMode) CALL DEBUG_CALL('CALC_GW', myThid )
725          IF (taveFreq.GT.0.) THEN  #endif
726            CALL TIMEAVE_CUMUL_1T(phiHydtave, phiHyd, Nr,           CALL TIMER_START('CALC_GW          [DYNAMICS]',myThid)
727       I                              deltaTclock, bi, bj, myThid)           CALL CALC_GW(
728            IF (ivdc_kappa.NE.0.) THEN       I                 bi,bj, KappaRU, KappaRV,
729              CALL TIMEAVE_CUMULATE(ConvectCountTave, ConvectCount, Nr,       I                 str13, str23, str33,
730       I                              deltaTclock, bi, bj, myThid)       I                 viscAh3d_00, viscAh3d_13, viscAh3d_23,
731            ENDIF       I                 myTime, myIter, myThid )
732          ENDIF          ENDIF
733  #endif /* ALLOW_TIMEAVE */          IF ( nonHydrostatic.OR.implicitIntGravWave )
734         &   CALL TIMESTEP_WVEL( bi,bj, myTime, myIter, myThid )
735            IF ( nonHydrostatic )
736         &   CALL TIMER_STOP ('CALC_GW          [DYNAMICS]',myThid)
737    #endif
738    
739    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
740    
741    C-    end of bi,bj loops
742         ENDDO         ENDDO
743        ENDDO        ENDDO
744    
745    #ifdef ALLOW_OBCS
746          IF (useOBCS) THEN
747            CALL OBCS_EXCHANGES( myThid )
748          ENDIF
749    #endif
750    
751    Cml(
752    C     In order to compare the variance of phiHydLow of a p/z-coordinate
753    C     run with etaH of a z/p-coordinate run the drift of phiHydLow
754    C     has to be removed by something like the following subroutine:
755    C      CALL REMOVE_MEAN_RL( 1, phiHydLow, maskInC, maskInC, rA, drF,
756    C     &                     'phiHydLow', myTime, myThid )
757    Cml)
758    
759    #ifdef ALLOW_DIAGNOSTICS
760          IF ( useDiagnostics ) THEN
761    
762           CALL DIAGNOSTICS_FILL(totPhihyd,'PHIHYD  ',0,Nr,0,1,1,myThid)
763           CALL DIAGNOSTICS_FILL(phiHydLow,'PHIBOT  ',0, 1,0,1,1,myThid)
764    
765           tmpFac = 1. _d 0
766           CALL DIAGNOSTICS_SCALE_FILL(totPhihyd,tmpFac,2,
767         &                                 'PHIHYDSQ',0,Nr,0,1,1,myThid)
768    
769           CALL DIAGNOSTICS_SCALE_FILL(phiHydLow,tmpFac,2,
770         &                                 'PHIBOTSQ',0, 1,0,1,1,myThid)
771    
772          ENDIF
773    #endif /* ALLOW_DIAGNOSTICS */
774    
775    #ifdef ALLOW_DEBUG
776          IF ( debugLevel .GE. debLevD ) THEN
777           CALL DEBUG_STATS_RL(1,EtaN,'EtaN (DYNAMICS)',myThid)
778           CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (DYNAMICS)',myThid)
779           CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (DYNAMICS)',myThid)
780           CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (DYNAMICS)',myThid)
781           CALL DEBUG_STATS_RL(Nr,theta,'Theta (DYNAMICS)',myThid)
782           CALL DEBUG_STATS_RL(Nr,salt,'Salt (DYNAMICS)',myThid)
783           CALL DEBUG_STATS_RL(Nr,gU,'Gu (DYNAMICS)',myThid)
784           CALL DEBUG_STATS_RL(Nr,gV,'Gv (DYNAMICS)',myThid)
785           CALL DEBUG_STATS_RL(Nr,gT,'Gt (DYNAMICS)',myThid)
786           CALL DEBUG_STATS_RL(Nr,gS,'Gs (DYNAMICS)',myThid)
787    #ifndef ALLOW_ADAMSBASHFORTH_3
788           CALL DEBUG_STATS_RL(Nr,guNm1,'GuNm1 (DYNAMICS)',myThid)
789           CALL DEBUG_STATS_RL(Nr,gvNm1,'GvNm1 (DYNAMICS)',myThid)
790           CALL DEBUG_STATS_RL(Nr,gtNm1,'GtNm1 (DYNAMICS)',myThid)
791           CALL DEBUG_STATS_RL(Nr,gsNm1,'GsNm1 (DYNAMICS)',myThid)
792    #endif
793          ENDIF
794    #endif
795    
796    #ifdef DYNAMICS_GUGV_EXCH_CHECK
797    C- jmc: For safety checking only: This Exchange here should not change
798    C       the solution. If solution changes, it means something is wrong,
799    C       but it does not mean that it is less wrong with this exchange.
800          IF ( debugLevel .GE. debLevE ) THEN
801           CALL EXCH_UV_XYZ_RL(gU,gV,.TRUE.,myThid)
802          ENDIF
803    #endif
804    
805    #ifdef ALLOW_DEBUG
806          IF (debugMode) CALL DEBUG_LEAVE( 'DYNAMICS', myThid )
807    #endif
808    
809        RETURN        RETURN
810        END        END
Legend:
Removed from v.1.68  
changed lines
  Added in v.1.167
  ViewVC Help
Powered by ViewVC 1.1.22