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