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

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

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

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