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