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

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

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

revision 1.74 by heimbach, Mon Jul 30 20:37:45 2001 UTC revision 1.122 by jmc, Sat Jul 30 23:39:48 2005 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    
4    #include "PACKAGES_CONFIG.h"
5  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
6    
7    CBOP
8    C     !ROUTINE: DYNAMICS
9    C     !INTERFACE:
10        SUBROUTINE DYNAMICS(myTime, myIter, myThid)        SUBROUTINE DYNAMICS(myTime, myIter, myThid)
11  C     /==========================================================\  C     !DESCRIPTION: \bv
12  C     | SUBROUTINE DYNAMICS                                      |  C     *==========================================================*
13  C     | o Controlling routine for the explicit part of the model |  C     | SUBROUTINE DYNAMICS                                      
14  C     |   dynamics.                                              |  C     | o Controlling routine for the explicit part of the model  
15  C     |==========================================================|  C     |   dynamics.                                              
16  C     | This routine evaluates the "dynamics" terms for each     |  C     *==========================================================*
17  C     | block of ocean in turn. Because the blocks of ocean have |  C     | This routine evaluates the "dynamics" terms for each      
18  C     | overlap regions they are independent of one another.     |  C     | block of ocean in turn. Because the blocks of ocean have  
19  C     | If terms involving lateral integrals are needed in this  |  C     | overlap regions they are independent of one another.      
20  C     | routine care will be needed. Similarly finite-difference |  C     | If terms involving lateral integrals are needed in this  
21  C     | operations with stencils wider than the overlap region   |  C     | routine care will be needed. Similarly finite-difference  
22  C     | require special consideration.                           |  C     | operations with stencils wider than the overlap region    
23  C     | Notes                                                    |  C     | require special consideration.                            
24  C     | =====                                                    |  C     | The algorithm...
25  C     | C*P* comments indicating place holders for which code is |  C     |
26  C     |      presently being developed.                          |  C     | "Correction Step"
27  C     \==========================================================/  C     | =================
28    C     | Here we update the horizontal velocities with the surface
29    C     | pressure such that the resulting flow is either consistent
30    C     | with the free-surface evolution or the rigid-lid:
31    C     |   U[n] = U* + dt x d/dx P
32    C     |   V[n] = V* + dt x d/dy P
33    C     |   W[n] = W* + dt x d/dz P  (NH mode)
34    C     |
35    C     | "Calculation of Gs"
36    C     | ===================
37    C     | This is where all the accelerations and tendencies (ie.
38    C     | physics, parameterizations etc...) are calculated
39    C     |   rho = rho ( theta[n], salt[n] )
40    C     |   b   = b(rho, theta)
41    C     |   K31 = K31 ( rho )
42    C     |   Gu[n] = Gu( u[n], v[n], wVel, b, ... )
43    C     |   Gv[n] = Gv( u[n], v[n], wVel, b, ... )
44    C     |   Gt[n] = Gt( theta[n], u[n], v[n], wVel, K31, ... )
45    C     |   Gs[n] = Gs( salt[n], u[n], v[n], wVel, K31, ... )
46    C     |
47    C     | "Time-stepping" or "Prediction"
48    C     | ================================
49    C     | The models variables are stepped forward with the appropriate
50    C     | time-stepping scheme (currently we use Adams-Bashforth II)
51    C     | - For momentum, the result is always *only* a "prediction"
52    C     | in that the flow may be divergent and will be "corrected"
53    C     | later with a surface pressure gradient.
54    C     | - Normally for tracers the result is the new field at time
55    C     | level [n+1} *BUT* in the case of implicit diffusion the result
56    C     | is also *only* a prediction.
57    C     | - We denote "predictors" with an asterisk (*).
58    C     |   U* = U[n] + dt x ( 3/2 Gu[n] - 1/2 Gu[n-1] )
59    C     |   V* = V[n] + dt x ( 3/2 Gv[n] - 1/2 Gv[n-1] )
60    C     |   theta[n+1] = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
61    C     |   salt[n+1] = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
62    C     | With implicit diffusion:
63    C     |   theta* = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
64    C     |   salt* = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
65    C     |   (1 + dt * K * d_zz) theta[n] = theta*
66    C     |   (1 + dt * K * d_zz) salt[n] = salt*
67    C     |
68    C     *==========================================================*
69    C     \ev
70    C     !USES:
71        IMPLICIT NONE        IMPLICIT NONE
   
72  C     == Global variables ===  C     == Global variables ===
73  #include "SIZE.h"  #include "SIZE.h"
74  #include "EEPARAMS.h"  #include "EEPARAMS.h"
75  #include "PARAMS.h"  #include "PARAMS.h"
76  #include "DYNVARS.h"  #include "DYNVARS.h"
77  #include "GRID.h"  #ifdef ALLOW_CD_CODE
78  #ifdef ALLOW_PASSIVE_TRACER  #include "CD_CODE_VARS.h"
 #include "TR1.h"  
79  #endif  #endif
80    #include "GRID.h"
81  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
82  # include "tamc.h"  # include "tamc.h"
83  # include "tamc_keys.h"  # include "tamc_keys.h"
84  # include "FFIELDS.h"  # include "FFIELDS.h"
85    # include "EOS.h"
86  # ifdef ALLOW_KPP  # ifdef ALLOW_KPP
87  #  include "KPP.h"  #  include "KPP.h"
88  # endif  # endif
 # ifdef ALLOW_GMREDI  
 #  include "GMREDI.h"  
 # endif  
89  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
90    
91  #ifdef ALLOW_TIMEAVE  C     !CALLING SEQUENCE:
92  #include "TIMEAVE_STATV.h"  C     DYNAMICS()
93  #endif  C      |
94    C      |-- CALC_EP_FORCING
95    C      |
96    C      |-- CALC_GRAD_PHI_SURF
97    C      |
98    C      |-- CALC_VISCOSITY
99    C      |
100    C      |-- CALC_PHI_HYD  
101    C      |
102    C      |-- MOM_FLUXFORM  
103    C      |
104    C      |-- MOM_VECINV    
105    C      |
106    C      |-- TIMESTEP      
107    C      |
108    C      |-- OBCS_APPLY_UV
109    C      |
110    C      |-- MOM_U_IMPLICIT_R      
111    C      |-- MOM_V_IMPLICIT_R      
112    C      |
113    C      |-- IMPLDIFF      
114    C      |
115    C      |-- OBCS_APPLY_UV
116    C      |
117    C      |-- CALC_GW
118    C      |
119    C      |-- DIAGNOSTICS_FILL
120    C      |-- DEBUG_STATS_RL
121    
122    C     !INPUT/OUTPUT PARAMETERS:
123  C     == Routine arguments ==  C     == Routine arguments ==
124  C     myTime - Current time in simulation  C     myTime - Current time in simulation
125  C     myIter - Current iteration number in simulation  C     myIter - Current iteration number in simulation
# Line 57  C     myThid - Thread number for this in Line 128  C     myThid - Thread number for this in
128        INTEGER myIter        INTEGER myIter
129        INTEGER myThid        INTEGER myThid
130    
131    C     !LOCAL VARIABLES:
132  C     == Local variables  C     == Local variables
133  C     xA, yA                 - Per block temporaries holding face areas  C     fVer[UV]               o fVer: Vertical flux term - note fVer
134  C     uTrans, vTrans, rTrans - Per block temporaries holding flow  C                                    is "pipelined" in the vertical
135  C                              transport  C                                    so we need an fVer for each
136  C                              o uTrans: Zonal transport  C                                    variable.
137  C                              o vTrans: Meridional transport  C     phiHydC    :: hydrostatic potential anomaly at cell center
138  C                              o rTrans: Vertical transport  C                   In z coords phiHyd is the hydrostatic potential
139  C     maskUp                   o maskUp: land/water mask for W points  C                      (=pressure/rho0) anomaly
140  C     fVer[STUV]               o fVer: Vertical flux term - note fVer  C                   In p coords phiHyd is the geopotential height anomaly.
141  C                                      is "pipelined" in the vertical  C     phiHydF    :: hydrostatic potential anomaly at middle between 2 centers
142  C                                      so we need an fVer for each  C     dPhiHydX,Y :: Gradient (X & Y directions) of hydrostatic potential anom.
143  C                                      variable.  C     phiSurfX,  ::  gradient of Surface potential (Pressure/rho, ocean)
144  C     rhoK, rhoKM1   - Density at current level, and level above  C     phiSurfY             or geopotential (atmos) in X and Y direction
145  C     phiHyd         - Hydrostatic part of the potential phiHydi.  C     guDissip   :: dissipation tendency (all explicit terms), u component
146  C                      In z coords phiHydiHyd is the hydrostatic  C     gvDissip   :: dissipation tendency (all explicit terms), v component
 C                      Potential (=pressure/rho0) anomaly  
 C                      In p coords phiHydiHyd is the geopotential  
 C                      surface height anomaly.  
 C     phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean)  
 C     phiSurfY             or geopotentiel (atmos) in X and Y direction  
 C     KappaRT,       - Total diffusion in vertical for T and S.  
 C     KappaRS          (background + spatially varying, isopycnal term).  
147  C     iMin, iMax     - Ranges and sub-block indices on which calculations  C     iMin, iMax     - Ranges and sub-block indices on which calculations
148  C     jMin, jMax       are applied.  C     jMin, jMax       are applied.
149  C     bi, bj  C     bi, bj
150  C     k, kup,        - Index for layer above and below. kup and kDown  C     k, kup,        - Index for layer above and below. kup and kDown
151  C     kDown, km1       are switched with layer to be the appropriate  C     kDown, km1       are switched with layer to be the appropriate
152  C                      index into fVerTerm.  C                      index into fVerTerm.
 C     tauAB - Adams-Bashforth timestepping weight: 0=forward ; 1/2=Adams-Bashf.  
       _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)  
       _RL fVerS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)  
       _RL fVerTr1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)  
153        _RL fVerU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
154        _RL fVerV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
155        _RL phiHyd  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL phiHydF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
156        _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiHydC (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
157        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL dPhiHydX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
158          _RL dPhiHydY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
159        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
160        _RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
161        _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL guDissip(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
162        _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL gvDissip(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
163        _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
164        _RL KappaRV (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRV (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
       _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
       _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
       _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
       _RL tauAB  
   
 C This is currently used by IVDC and Diagnostics  
       _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
165    
166        INTEGER iMin, iMax        INTEGER iMin, iMax
167        INTEGER jMin, jMax        INTEGER jMin, jMax
168        INTEGER bi, bj        INTEGER bi, bj
169        INTEGER i, j        INTEGER i, j
170        INTEGER k, km1, kup, kDown        INTEGER k, km1, kp1, kup, kDown
171    
172    #ifdef ALLOW_DIAGNOSTICS
173          _RL tmpFac
174    #endif /* ALLOW_DIAGNOSTICS */
175    
 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)  
176    
177  C---    The algorithm...  C---    The algorithm...
178  C  C
# Line 169  C         salt* = salt[n] + dt x ( 3/2 G Line 217  C         salt* = salt[n] + dt x ( 3/2 G
217  C         (1 + dt * K * d_zz) theta[n] = theta*  C         (1 + dt * K * d_zz) theta[n] = theta*
218  C         (1 + dt * K * d_zz) salt[n] = salt*  C         (1 + dt * K * d_zz) salt[n] = salt*
219  C---  C---
220    CEOP
221    
222  #ifdef ALLOW_AUTODIFF_TAMC  C-- Call to routine for calculation of
223  C--   dummy statement to end declaration part  C   Eliassen-Palm-flux-forced U-tendency,
224        ikey = 1  C   if desired:
225  #endif /* ALLOW_AUTODIFF_TAMC */  #ifdef INCLUDE_EP_FORCING_CODE
226          CALL CALC_EP_FORCING(myThid)
227  C--   Set up work arrays with valid (i.e. not NaN) values  #endif
 C     These inital values do not alter the numerical results. They  
 C     just ensure that all memory references are to valid floating  
 C     point numbers. This prevents spurious hardware signals due to  
 C     uninitialised but inert locations.  
       DO j=1-OLy,sNy+OLy  
        DO i=1-OLx,sNx+OLx  
         xA(i,j)      = 0. _d 0  
         yA(i,j)      = 0. _d 0  
         uTrans(i,j)  = 0. _d 0  
         vTrans(i,j)  = 0. _d 0  
         DO k=1,Nr  
          phiHyd(i,j,k)  = 0. _d 0  
          KappaRU(i,j,k) = 0. _d 0  
          KappaRV(i,j,k) = 0. _d 0  
          sigmaX(i,j,k) = 0. _d 0  
          sigmaY(i,j,k) = 0. _d 0  
          sigmaR(i,j,k) = 0. _d 0  
         ENDDO  
         rhoKM1 (i,j) = 0. _d 0  
         rhok   (i,j) = 0. _d 0  
         phiSurfX(i,j) = 0. _d 0  
         phiSurfY(i,j) = 0. _d 0  
        ENDDO  
       ENDDO  
   
228    
229  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
230  C--   HPF directive to help TAMC  C--   HPF directive to help TAMC
# Line 211  CHPF$ INDEPENDENT Line 235  CHPF$ INDEPENDENT
235    
236  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
237  C--    HPF directive to help TAMC  C--    HPF directive to help TAMC
238  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS,fVerU,fVerV  CHPF$  INDEPENDENT, NEW (fVerU,fVerV
239  CHPF$&                  ,phiHyd,utrans,vtrans,xA,yA  CHPF$&                  ,phiHydF
240  CHPF$&                  ,KappaRT,KappaRS,KappaRU,KappaRV  CHPF$&                  ,KappaRU,KappaRV
241  CHPF$&                  )  CHPF$&                  )
242  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
243    
# Line 222  CHPF$&                  ) Line 246  CHPF$&                  )
246  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
247            act1 = bi - myBxLo(myThid)            act1 = bi - myBxLo(myThid)
248            max1 = myBxHi(myThid) - myBxLo(myThid) + 1            max1 = myBxHi(myThid) - myBxLo(myThid) + 1
   
249            act2 = bj - myByLo(myThid)            act2 = bj - myByLo(myThid)
250            max2 = myByHi(myThid) - myByLo(myThid) + 1            max2 = myByHi(myThid) - myByLo(myThid) + 1
   
251            act3 = myThid - 1            act3 = myThid - 1
252            max3 = nTx*nTy            max3 = nTx*nTy
   
253            act4 = ikey_dynamics - 1            act4 = ikey_dynamics - 1
254              idynkey = (act1 + 1) + act2*max1
           ikey = (act1 + 1) + act2*max1  
255       &                      + act3*max1*max2       &                      + act3*max1*max2
256       &                      + act4*max1*max2*max3       &                      + act4*max1*max2*max3
257  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
258    
259  C--     Set up work arrays that need valid initial values  C--   Set up work arrays with valid (i.e. not NaN) values
260          DO j=1-OLy,sNy+OLy  C     These inital values do not alter the numerical results. They
261           DO i=1-OLx,sNx+OLx  C     just ensure that all memory references are to valid floating
262            rTrans (i,j)   = 0. _d 0  C     point numbers. This prevents spurious hardware signals due to
263            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  
           fVerTr1(i,j,1) = 0. _d 0  
           fVerTr1(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  
264    
265          DO k=1,Nr          DO k=1,Nr
266           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
267            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
268  C This is currently also used by IVDC and Diagnostics             KappaRU(i,j,k) = 0. _d 0
269             ConvectCount(i,j,k) = 0.             KappaRV(i,j,k) = 0. _d 0
270             KappaRT(i,j,k) = 0. _d 0  #ifdef ALLOW_AUTODIFF_TAMC
271             KappaRS(i,j,k) = 0. _d 0  cph(
272    c--   need some re-initialisation here to break dependencies
273    cph)
274               gU(i,j,k,bi,bj) = 0. _d 0
275               gV(i,j,k,bi,bj) = 0. _d 0
276    #endif
277            ENDDO            ENDDO
278           ENDDO           ENDDO
279          ENDDO          ENDDO
280            DO j=1-OLy,sNy+OLy
281          iMin = 1-OLx+1           DO i=1-OLx,sNx+OLx
282          iMax = sNx+OLx            fVerU  (i,j,1) = 0. _d 0
283          jMin = 1-OLy+1            fVerU  (i,j,2) = 0. _d 0
284          jMax = sNy+OLy            fVerV  (i,j,1) = 0. _d 0
285              fVerV  (i,j,2) = 0. _d 0
286              phiHydF (i,j)  = 0. _d 0
287  #ifdef ALLOW_AUTODIFF_TAMC            phiHydC (i,j)  = 0. _d 0
288  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte            dPhiHydX(i,j)  = 0. _d 0
289  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte            dPhiHydY(i,j)  = 0. _d 0
290  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte            phiSurfX(i,j)  = 0. _d 0
291  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte            phiSurfY(i,j)  = 0. _d 0
292  #ifdef ALLOW_PASSIVE_TRACER            guDissip(i,j)  = 0. _d 0
293  CADJ STORE tr1  (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte            gvDissip(i,j)  = 0. _d 0
294  #endif           ENDDO
 #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)  
295          ENDDO          ENDDO
296    
297  #ifdef ALLOW_AUTODIFF_TAMC  C--     Start computation of dynamics
298  cph avoids recomputation of integrate_for_w          iMin = 0
299  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte          iMax = sNx+1
300  #endif /* ALLOW_AUTODIFF_TAMC */          jMin = 0
301            jMax = sNy+1
 #ifdef  ALLOW_OBCS  
 C--     Calculate future values on open boundaries  
         IF (useOBCS) THEN  
           CALL OBCS_CALC( bi, bj, myTime+deltaT,  
      I            uVel, vVel, wVel, theta, salt,  
      I            myThid )  
         ENDIF  
 #endif  /* ALLOW_OBCS */  
   
 C--     Determines forcing terms based on external fields  
 C       relaxation terms, etc.  
         CALL EXTERNAL_FORCING_SURF(  
      I             bi, bj, iMin, iMax, jMin, jMax,  
      I             myThid )  
 #ifdef ALLOW_AUTODIFF_TAMC  
 cph needed for KPP  
 CADJ STORE surfacetendencyU(:,:,bi,bj)  
 CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE surfacetendencyV(:,:,bi,bj)  
 CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE surfacetendencyS(:,:,bi,bj)  
 CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE surfacetendencyT(:,:,bi,bj)  
 CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 #ifdef  ALLOW_GMREDI  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE sigmaX(:,:,:) = comlev1, key=ikey, byte=isbyte  
 CADJ STORE sigmaY(:,:,:) = comlev1, key=ikey, byte=isbyte  
 CADJ STORE sigmaR(:,:,:) = comlev1, key=ikey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
 C--     Calculate iso-neutral slopes for the GM/Redi parameterisation  
         IF (useGMRedi) THEN  
           DO k=1,Nr  
             CALL GMREDI_CALC_TENSOR(  
      I             bi, bj, iMin, iMax, jMin, jMax, k,  
      I             sigmaX, sigmaY, sigmaR,  
      I             myThid )  
           ENDDO  
 #ifdef ALLOW_AUTODIFF_TAMC  
         ELSE  
           DO k=1, Nr  
             CALL GMREDI_CALC_TENSOR_DUMMY(  
      I             bi, bj, iMin, iMax, jMin, jMax, k,  
      I             sigmaX, sigmaY, sigmaR,  
      I             myThid )  
           ENDDO  
 #endif /* ALLOW_AUTODIFF_TAMC */  
         ENDIF  
302    
303  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
304  CADJ STORE Kwx(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE wvel (:,:,:,bi,bj) =
305  CADJ STORE Kwy(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte  CADJ &     comlev1_bibj, key = idynkey, byte = isbyte
 CADJ STORE Kwz(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte  
306  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
307    
308  #endif  /* ALLOW_GMREDI */  C--     Explicit part of the Surface Potentiel Gradient (add in TIMESTEP)
309    C       (note: this loop will be replaced by CALL CALC_GRAD_ETA)
310  #ifdef  ALLOW_KPP          IF (implicSurfPress.NE.1.) THEN
311  C--     Compute KPP mixing coefficients            CALL CALC_GRAD_PHI_SURF(
312          IF (useKPP) THEN       I         bi,bj,iMin,iMax,jMin,jMax,
313            CALL KPP_CALC(       I         etaN,
314       I                  bi, bj, myTime, myThid )       O         phiSurfX,phiSurfY,
315  #ifdef ALLOW_AUTODIFF_TAMC       I         myThid )                        
         ELSE  
           CALL KPP_CALC_DUMMY(  
      I                  bi, bj, myTime, myThid )  
 #endif /* ALLOW_AUTODIFF_TAMC */  
316          ENDIF          ENDIF
317    
318  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
319  CADJ STORE KPPghat   (:,:,:,bi,bj)  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=idynkey, byte=isbyte
320  CADJ &   , KPPviscAz (:,:,:,bi,bj)  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=idynkey, byte=isbyte
321  CADJ &   , KPPdiffKzT(:,:,:,bi,bj)  #ifdef ALLOW_KPP
322  CADJ &   , KPPdiffKzS(:,:,:,bi,bj)  CADJ STORE KPPviscAz (:,:,:,bi,bj)
323  CADJ &   , KPPfrac   (:,:  ,bi,bj)  CADJ &                 = comlev1_bibj, key=idynkey, byte=isbyte
324  CADJ &                 = comlev1_bibj, key=ikey, byte=isbyte  #endif /* ALLOW_KPP */
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 #endif  /* ALLOW_KPP */  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE KappaRT(:,:,:)     = comlev1_bibj, key = ikey, byte = isbyte  
 CADJ STORE KappaRS(:,:,:)     = comlev1_bibj, key = ikey, byte = isbyte  
 CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
 CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
 CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
 CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
 #ifdef ALLOW_PASSIVE_TRACER  
 CADJ STORE tr1  (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
 #endif  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 #ifdef ALLOW_AIM  
 C       AIM - atmospheric intermediate model, physics package code.  
 C note(jmc) : phiHyd=0 at this point but is not really used in Molteni Physics  
         IF ( useAIM ) THEN  
          CALL TIMER_START('AIM_DO_ATMOS_PHYS      [DYNAMICS]', myThid)  
          CALL AIM_DO_ATMOS_PHYSICS( phiHyd, bi, bj, myTime, myThid )  
          CALL TIMER_STOP ('AIM_DO_ATMOS_PHYS      [DYNAMICS]', myThid)  
         ENDIF  
 #endif /* ALLOW_AIM */  
   
   
 C--     Start of thermodynamics loop  
         DO k=Nr,1,-1  
 #ifdef ALLOW_AUTODIFF_TAMC  
 C? Patrick Is this formula correct?  
 cph Yes, but I rewrote it.  
 cph Also, the KappaR? need the index and subscript k!  
          kkey = (ikey-1)*Nr + k  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 C--       km1    Points to level above k (=k-1)  
 C--       kup    Cycles through 1,2 to point to layer above  
 C--       kDown  Cycles through 2,1 to point to current layer  
   
           km1  = MAX(1,k-1)  
           kup  = 1+MOD(k+1,2)  
           kDown= 1+MOD(k,2)  
   
           iMin = 1-OLx  
           iMax = sNx+OLx  
           jMin = 1-OLy  
           jMax = sNy+OLy  
   
 C--       Get temporary terms used by tendency routines  
           CALL CALC_COMMON_FACTORS (  
      I         bi,bj,iMin,iMax,jMin,jMax,k,  
      O         xA,yA,uTrans,vTrans,rTrans,maskUp,  
      I         myThid)  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE KappaRT(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte  
 CADJ STORE KappaRS(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte  
325  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
326    
327  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL
328  C--      Calculate the total vertical diffusivity  C--      Calculate the total vertical diffusivity
329           CALL CALC_DIFFUSIVITY(          DO k=1,Nr
330             CALL CALC_VISCOSITY(
331       I        bi,bj,iMin,iMax,jMin,jMax,k,       I        bi,bj,iMin,iMax,jMin,jMax,k,
332       I        maskUp,       O        KappaRU,KappaRV,
      O        KappaRT,KappaRS,KappaRU,KappaRV,  
333       I        myThid)       I        myThid)
334           ENDDO
335  #endif  #endif
336    
           iMin = 1-OLx+2  
           iMax = sNx+OLx-1  
           jMin = 1-OLy+2  
           jMax = sNy+OLy-1  
   
 C--      Calculate active tracer tendencies (gT,gS,...)  
 C        and step forward storing result in gTnm1, gSnm1, etc.  
          IF ( tempStepping ) THEN  
            CALL CALC_GT(  
      I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,  
      I         xA,yA,uTrans,vTrans,rTrans,maskUp,  
      I         KappaRT,  
      U         fVerT,  
      I         myTime, myThid)  
            tauAB = 0.5d0 + abEps  
            CALL TIMESTEP_TRACER(  
      I         bi,bj,iMin,iMax,jMin,jMax,k,tauAB,  
      I         theta, gT,  
      U         gTnm1,  
      I         myIter, myThid)  
          ENDIF  
          IF ( saltStepping ) THEN  
            CALL CALC_GS(  
      I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,  
      I         xA,yA,uTrans,vTrans,rTrans,maskUp,  
      I         KappaRS,  
      U         fVerS,  
      I         myTime, myThid)  
            tauAB = 0.5d0 + abEps  
            CALL TIMESTEP_TRACER(  
      I         bi,bj,iMin,iMax,jMin,jMax,k,tauAB,  
      I         salt, gS,  
      U         gSnm1,  
      I         myIter, myThid)  
          ENDIF  
 #ifdef ALLOW_PASSIVE_TRACER  
          IF ( tr1Stepping ) THEN  
            CALL CALC_GTR1(  
      I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,  
      I         xA,yA,uTrans,vTrans,rTrans,maskUp,  
      I         KappaRT,  
      U         fVerTr1,  
      I         myTime, myThid)  
            tauAB = 0.5d0 + abEps  
            CALL TIMESTEP_TRACER(  
      I         bi,bj,iMin,iMax,jMin,jMax,k,tauAB,  
      I         Tr1, gTr1,  
      U         gTr1NM1,  
      I         myIter, myThid)  
          ENDIF  
 #endif  
   
 #ifdef   ALLOW_OBCS  
 C--      Apply open boundary conditions  
          IF (useOBCS) THEN  
            CALL OBCS_APPLY_TS( bi, bj, k, gTnm1, gSnm1, myThid )  
          END IF  
 #endif   /* ALLOW_OBCS */  
   
 C--      Freeze water  
          IF (allowFreezing) THEN  
337  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
338  CADJ STORE gTNm1(:,:,k,bi,bj) = comlev1_bibj_k  CADJ STORE KappaRU(:,:,:)
339  CADJ &   , key = kkey, byte = isbyte  CADJ &                 = comlev1_bibj, key=idynkey, byte=isbyte
340    CADJ STORE KappaRV(:,:,:)
341    CADJ &                 = comlev1_bibj, key=idynkey, byte=isbyte
342  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
             CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )  
          END IF  
   
 C--     end of thermodynamic k loop (Nr:1)  
         ENDDO  
   
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 C? Patrick? What about this one?  
 cph Keys iikey and idkey don't seem to be needed  
 cph since storing occurs on different tape for each  
 cph impldiff call anyways.  
 cph Thus, common block comlev1_impl isn't needed either.  
 cph Storing below needed in the case useGMREDI.  
         iikey = (ikey-1)*maximpl  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 C--     Implicit diffusion  
         IF (implicitDiffusion) THEN  
   
          IF (tempStepping) THEN  
 #ifdef ALLOW_AUTODIFF_TAMC  
             idkey = iikey + 1  
 CADJ STORE gTNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
             CALL IMPLDIFF(  
      I         bi, bj, iMin, iMax, jMin, jMax,  
      I         deltaTtracer, KappaRT, recip_HFacC,  
      U         gTNm1,  
      I         myThid )  
          ENDIF  
   
          IF (saltStepping) THEN  
 #ifdef ALLOW_AUTODIFF_TAMC  
          idkey = iikey + 2  
 CADJ STORE gSNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
             CALL IMPLDIFF(  
      I         bi, bj, iMin, iMax, jMin, jMax,  
      I         deltaTtracer, KappaRS, recip_HFacC,  
      U         gSNm1,  
      I         myThid )  
          ENDIF  
   
 #ifdef ALLOW_PASSIVE_TRACER  
          IF (tr1Stepping) THEN  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE gTr1Nm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
           CALL IMPLDIFF(  
      I      bi, bj, iMin, iMax, jMin, jMax,  
      I      deltaTtracer, KappaRT, recip_HFacC,  
      U      gTr1Nm1,  
      I      myThid )  
          ENDIF  
 #endif  
   
 #ifdef   ALLOW_OBCS  
 C--      Apply open boundary conditions  
          IF (useOBCS) THEN  
            DO K=1,Nr  
              CALL OBCS_APPLY_TS( bi, bj, k, gTnm1, gSnm1, myThid )  
            ENDDO  
          END IF  
 #endif   /* ALLOW_OBCS */  
   
 C--     End If implicitDiffusion  
         ENDIF  
   
 C--     Start computation of dynamics  
         iMin = 1-OLx+2  
         iMax = sNx+OLx-1  
         jMin = 1-OLy+2  
         jMax = sNy+OLy-1  
   
 C--     Explicit part of the Surface Potentiel Gradient (add in TIMESTEP)  
 C       (note: this loop will be replaced by CALL CALC_GRAD_ETA)  
         IF (implicSurfPress.NE.1.) THEN  
           CALL CALC_GRAD_PHI_SURF(  
      I         bi,bj,iMin,iMax,jMin,jMax,  
      I         etaN,  
      O         phiSurfX,phiSurfY,  
      I         myThid )                          
         ENDIF  
343    
344  C--     Start of dynamics loop  C--     Start of dynamics loop
345          DO k=1,Nr          DO k=1,Nr
# Line 663  C--       kup    Cycles through 1,2 to p Line 349  C--       kup    Cycles through 1,2 to p
349  C--       kDown  Cycles through 2,1 to point to current layer  C--       kDown  Cycles through 2,1 to point to current layer
350    
351            km1  = MAX(1,k-1)            km1  = MAX(1,k-1)
352              kp1  = MIN(k+1,Nr)
353            kup  = 1+MOD(k+1,2)            kup  = 1+MOD(k+1,2)
354            kDown= 1+MOD(k,2)            kDown= 1+MOD(k,2)
355    
356    #ifdef ALLOW_AUTODIFF_TAMC
357             kkey = (idynkey-1)*Nr + k
358    c
359    CADJ STORE totphihyd (:,:,k,bi,bj)
360    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
361    CADJ STORE theta (:,:,k,bi,bj)
362    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
363    CADJ STORE salt  (:,:,k,bi,bj)
364    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
365    #endif /* ALLOW_AUTODIFF_TAMC */
366    
367  C--      Integrate hydrostatic balance for phiHyd with BC of  C--      Integrate hydrostatic balance for phiHyd with BC of
368  C        phiHyd(z=0)=0  C        phiHyd(z=0)=0
369  C        distinguishe between Stagger and Non Stagger time stepping           CALL CALC_PHI_HYD(
          IF (staggerTimeStep) THEN  
            CALL CALC_PHI_HYD(  
      I        bi,bj,iMin,iMax,jMin,jMax,k,  
      I        gTnm1, gSnm1,  
      U        phiHyd,  
      I        myThid )  
          ELSE  
            CALL CALC_PHI_HYD(  
370       I        bi,bj,iMin,iMax,jMin,jMax,k,       I        bi,bj,iMin,iMax,jMin,jMax,k,
371       I        theta, salt,       I        theta, salt,
372       U        phiHyd,       U        phiHydF,
373       I        myThid )       O        phiHydC, dPhiHydX, dPhiHydY,
374           ENDIF       I        myTime, myIter, myThid )
375    
376  C--      Calculate accelerations in the momentum equations (gU, gV, ...)  C--      Calculate accelerations in the momentum equations (gU, gV, ...)
377  C        and step forward storing the result in gUnm1, gVnm1, etc...  C        and step forward storing the result in gU, gV, etc...
378           IF ( momStepping ) THEN           IF ( momStepping ) THEN
379             CALL CALC_MOM_RHS(  #ifdef ALLOW_MOM_FLUXFORM
380               IF (.NOT. vectorInvariantMomentum) CALL MOM_FLUXFORM(
381       I         bi,bj,iMin,iMax,jMin,jMax,k,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax,k,kup,kDown,
382       I         phiHyd,KappaRU,KappaRV,       I         KappaRU, KappaRV,
383       U         fVerU, fVerV,       U         fVerU, fVerV,
384       I         myTime, myThid)       O         guDissip, gvDissip,
385         I         myTime, myIter, myThid)
386    #endif
387    #ifdef ALLOW_MOM_VECINV
388               IF (vectorInvariantMomentum) CALL MOM_VECINV(
389         I         bi,bj,iMin,iMax,jMin,jMax,k,kup,kDown,
390         I         KappaRU, KappaRV,
391         U         fVerU, fVerV,
392         O         guDissip, gvDissip,
393         I         myTime, myIter, myThid)
394    #endif
395             CALL TIMESTEP(             CALL TIMESTEP(
396       I         bi,bj,iMin,iMax,jMin,jMax,k,       I         bi,bj,iMin,iMax,jMin,jMax,k,
397       I         phiHyd, phiSurfX, phiSurfY,       I         dPhiHydX,dPhiHydY, phiSurfX, phiSurfY,
398       I         myIter, myThid)       I         guDissip, gvDissip,
399         I         myTime, myIter, myThid)
400    
401  #ifdef   ALLOW_OBCS  #ifdef   ALLOW_OBCS
402  C--      Apply open boundary conditions  C--      Apply open boundary conditions
403           IF (useOBCS) THEN             IF (useOBCS) THEN
404             CALL OBCS_APPLY_UV( bi, bj, k, gUnm1, gVnm1, myThid )               CALL OBCS_APPLY_UV( bi, bj, k, gU, gV, myThid )
405           END IF             ENDIF
406  #endif   /* ALLOW_OBCS */  #endif   /* ALLOW_OBCS */
407    
 #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 */  
408           ENDIF           ENDIF
409    
410    
411  C--     end of dynamics k loop (1:Nr)  C--     end of dynamics k loop (1:Nr)
412          ENDDO          ENDDO
413    
414    C--     Implicit Vertical advection & viscosity
415    #ifdef INCLUDE_IMPLVERTADV_CODE
416  C--     Implicit viscosity          IF ( momImplVertAdv ) THEN
417          IF (implicitViscosity.AND.momStepping) THEN            CALL MOM_U_IMPLICIT_R( kappaRU,
418         I                           bi, bj, myTime, myIter, myThid )
419              CALL MOM_V_IMPLICIT_R( kappaRV,
420         I                           bi, bj, myTime, myIter, myThid )
421            ELSEIF ( implicitViscosity ) THEN
422    #else /* INCLUDE_IMPLVERTADV_CODE */
423            IF     ( implicitViscosity ) THEN
424    #endif /* INCLUDE_IMPLVERTADV_CODE */
425  #ifdef    ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
426            idkey = iikey + 3  CADJ STORE KappaRU(:,:,:) = comlev1_bibj , key=idynkey, byte=isbyte
427  CADJ STORE gUNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE gU(:,:,:,bi,bj) = comlev1_bibj , key=idynkey, byte=isbyte
428  #endif    /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
429            CALL IMPLDIFF(            CALL IMPLDIFF(
430       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
431       I         deltaTmom, KappaRU,recip_HFacW,       I         0, KappaRU,recip_HFacW,
432       U         gUNm1,       U         gU,
433       I         myThid )       I         myThid )
434  #ifdef    ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
435            idkey = iikey + 4  CADJ STORE KappaRV(:,:,:) = comlev1_bibj , key=idynkey, byte=isbyte
436  CADJ STORE gVNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE gV(:,:,:,bi,bj) = comlev1_bibj , key=idynkey, byte=isbyte
437  #endif    /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
438            CALL IMPLDIFF(            CALL IMPLDIFF(
439       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
440       I         deltaTmom, KappaRV,recip_HFacS,       I         0, KappaRV,recip_HFacS,
441       U         gVNm1,       U         gV,
442       I         myThid )       I         myThid )
443            ENDIF
444    
445  #ifdef   ALLOW_OBCS  #ifdef   ALLOW_OBCS
446  C--      Apply open boundary conditions  C--      Apply open boundary conditions
447           IF (useOBCS) THEN          IF ( useOBCS .AND.(implicitViscosity.OR.momImplVertAdv) ) THEN
448             DO K=1,Nr             DO K=1,Nr
449               CALL OBCS_APPLY_UV( bi, bj, k, gUnm1, gVnm1, myThid )               CALL OBCS_APPLY_UV( bi, bj, k, gU, gV, myThid )
450             ENDDO             ENDDO
451           END IF          ENDIF
452  #endif   /* ALLOW_OBCS */  #endif   /* ALLOW_OBCS */
453    
454  #ifdef    INCLUDE_CD_CODE  #ifdef    ALLOW_CD_CODE
455            IF (implicitViscosity.AND.useCDscheme) THEN
456  #ifdef    ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
457            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  
458  #endif    /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
459            CALL IMPLDIFF(            CALL IMPLDIFF(
460       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
461       I         deltaTmom, KappaRU,recip_HFacW,       I         0, KappaRU,recip_HFacW,
462       U         vVelD,       U         vVelD,
463       I         myThid )       I         myThid )
464  #ifdef    ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
465            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  
466  #endif    /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
467            CALL IMPLDIFF(            CALL IMPLDIFF(
468       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
469       I         deltaTmom, KappaRV,recip_HFacS,       I         0, KappaRV,recip_HFacS,
470       U         uVelD,       U         uVelD,
471       I         myThid )       I         myThid )
 #endif    /* INCLUDE_CD_CODE */  
 C--     End If implicitViscosity.AND.momStepping  
472          ENDIF          ENDIF
473    #endif    /* ALLOW_CD_CODE */
474    C--     End implicit Vertical advection & viscosity
475    
 Cjmc : add for phiHyd output <- but not working if multi tile per CPU  
 c       IF ( DIFFERENT_MULTIPLE(dumpFreq,myTime+deltaTClock,myTime)  
 c    &  .AND. buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN  
 c         WRITE(suff,'(I10.10)') myIter+1  
 c         CALL WRITE_FLD_XYZ_RL('PH.',suff,phiHyd,myIter+1,myThid)  
 c       ENDIF  
 Cjmc(end)  
   
 #ifdef ALLOW_TIMEAVE  
         IF (taveFreq.GT.0.) THEN  
           CALL TIMEAVE_CUMUL_1T(phiHydtave, phiHyd, Nr,  
      I                              deltaTclock, bi, bj, myThid)  
           IF (ivdc_kappa.NE.0.) THEN  
             CALL TIMEAVE_CUMULATE(ConvectCountTave, ConvectCount, Nr,  
      I                              deltaTclock, bi, bj, myThid)  
           ENDIF  
         ENDIF  
 #endif /* ALLOW_TIMEAVE */  
   
476         ENDDO         ENDDO
477        ENDDO        ENDDO
478    
479  #ifndef EXCLUDE_DEBUGMODE  #ifdef ALLOW_OBCS
480        If (debugMode) THEN        IF (useOBCS) THEN
481           CALL OBCS_PRESCRIBE_EXCHANGES(myThid)
482          ENDIF
483    #endif
484    
485    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
486    
487    #ifdef ALLOW_NONHYDROSTATIC
488    C--   Step forward W field in N-H algorithm
489          IF ( momStepping .AND. nonHydrostatic ) THEN
490    #ifdef ALLOW_DEBUG
491              IF ( debugLevel .GE. debLevB )
492         &     CALL DEBUG_CALL('CALC_GW',myThid)
493    #endif
494             CALL TIMER_START('CALC_GW          [DYNAMICS]',myThid)
495             CALL CALC_GW( myTime, myIter, myThid )
496             CALL TIMER_STOP ('CALC_GW          [DYNAMICS]',myThid)
497          ENDIF
498    #endif
499    
500    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
501    
502    Cml(
503    C     In order to compare the variance of phiHydLow of a p/z-coordinate
504    C     run with etaH of a z/p-coordinate run the drift of phiHydLow
505    C     has to be removed by something like the following subroutine:
506    C      CALL REMOVE_MEAN_RL( 1, phiHydLow, maskH, maskH, rA, drF,
507    C     &                'phiHydLow', myThid )
508    Cml)
509    
510    #ifdef ALLOW_DIAGNOSTICS
511          IF ( usediagnostics ) THEN
512    
513           CALL DIAGNOSTICS_FILL(totPhihyd,'PHIHYD  ',0,Nr,0,1,1,myThid)
514           CALL DIAGNOSTICS_FILL(phiHydLow,'PHIBOT  ',0, 1,0,1,1,myThid)
515    
516           tmpFac = 1. _d 0
517           CALL DIAGNOSTICS_SCALE_FILL(totPhihyd,tmpFac,2,
518         &                                 'PHIHYDSQ',0,Nr,0,1,1,myThid)
519    
520           CALL DIAGNOSTICS_SCALE_FILL(phiHydLow,tmpFac,2,
521         &                                 'PHIBOTSQ',0, 1,0,1,1,myThid)
522    
523          ENDIF
524    #endif /* ALLOW_DIAGNOSTICS */
525          
526    #ifdef ALLOW_DEBUG
527          If ( debugLevel .GE. debLevB ) THEN
528         CALL DEBUG_STATS_RL(1,EtaN,'EtaN (DYNAMICS)',myThid)         CALL DEBUG_STATS_RL(1,EtaN,'EtaN (DYNAMICS)',myThid)
529         CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (DYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (DYNAMICS)',myThid)
530         CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (DYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (DYNAMICS)',myThid)
531         CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (DYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (DYNAMICS)',myThid)
532         CALL DEBUG_STATS_RL(Nr,theta,'Theta (DYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,theta,'Theta (DYNAMICS)',myThid)
533         CALL DEBUG_STATS_RL(Nr,salt,'Salt (DYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,salt,'Salt (DYNAMICS)',myThid)
534         CALL DEBUG_STATS_RL(Nr,Gu,'Gu (DYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,gU,'Gu (DYNAMICS)',myThid)
535         CALL DEBUG_STATS_RL(Nr,Gv,'Gv (DYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,gV,'Gv (DYNAMICS)',myThid)
536         CALL DEBUG_STATS_RL(Nr,Gt,'Gt (DYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,gT,'Gt (DYNAMICS)',myThid)
537         CALL DEBUG_STATS_RL(Nr,Gs,'Gs (DYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,gS,'Gs (DYNAMICS)',myThid)
538         CALL DEBUG_STATS_RL(Nr,GuNm1,'GuNm1 (DYNAMICS)',myThid)  #ifndef ALLOW_ADAMSBASHFORTH_3
539         CALL DEBUG_STATS_RL(Nr,GvNm1,'GvNm1 (DYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,guNm1,'GuNm1 (DYNAMICS)',myThid)
540         CALL DEBUG_STATS_RL(Nr,GtNm1,'GtNm1 (DYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,gvNm1,'GvNm1 (DYNAMICS)',myThid)
541         CALL DEBUG_STATS_RL(Nr,GsNm1,'GsNm1 (DYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,gtNm1,'GtNm1 (DYNAMICS)',myThid)
542           CALL DEBUG_STATS_RL(Nr,gsNm1,'GsNm1 (DYNAMICS)',myThid)
543    #endif
544        ENDIF        ENDIF
545  #endif  #endif
546    
Legend:
Removed from v.1.74  
changed lines
  Added in v.1.122
  ViewVC Help
Powered by ViewVC 1.1.22