/[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.66 by heimbach, Sun Mar 25 22:33:52 2001 UTC revision 1.131 by heimbach, Wed Mar 29 17:00:39 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  #endif /* ALLOW_AUTODIFF_TAMC */  # 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_PTRACERS
95  #ifdef ALLOW_TIMEAVE  #  include "PTRACERS_SIZE.h"
96  #include "TIMEAVE_STATV.h"  #  include "PTRACERS.h"
97  #endif  # endif
98    # ifdef ALLOW_OBCS
99    #  include "OBCS.h"
100    #  ifdef ALLOW_PTRACERS
101    #   include "OBCS_PTRACERS.h"
102    #  endif
103    # endif
104    #endif /* ALLOW_AUTODIFF_TAMC */
105    
106    C     !CALLING SEQUENCE:
107    C     DYNAMICS()
108    C      |
109    C      |-- CALC_EP_FORCING
110    C      |
111    C      |-- CALC_GRAD_PHI_SURF
112    C      |
113    C      |-- CALC_VISCOSITY
114    C      |
115    C      |-- CALC_PHI_HYD  
116    C      |
117    C      |-- MOM_FLUXFORM  
118    C      |
119    C      |-- MOM_VECINV    
120    C      |
121    C      |-- TIMESTEP      
122    C      |
123    C      |-- OBCS_APPLY_UV
124    C      |
125    C      |-- MOM_U_IMPLICIT_R      
126    C      |-- MOM_V_IMPLICIT_R      
127    C      |
128    C      |-- IMPLDIFF      
129    C      |
130    C      |-- OBCS_APPLY_UV
131    C      |
132    C      |-- CALC_GW
133    C      |
134    C      |-- DIAGNOSTICS_FILL
135    C      |-- DEBUG_STATS_RL
136    
137    C     !INPUT/OUTPUT PARAMETERS:
138  C     == Routine arguments ==  C     == Routine arguments ==
139  C     myTime - Current time in simulation  C     myTime - Current time in simulation
140  C     myIter - Current iteration number in simulation  C     myIter - Current iteration number in simulation
# Line 51  C     myThid - Thread number for this in Line 143  C     myThid - Thread number for this in
143        INTEGER myIter        INTEGER myIter
144        INTEGER myThid        INTEGER myThid
145    
146    C     !LOCAL VARIABLES:
147  C     == Local variables  C     == Local variables
148  C     xA, yA                 - Per block temporaries holding face areas  C     fVer[UV]               o fVer: Vertical flux term - note fVer
149  C     uTrans, vTrans, rTrans - Per block temporaries holding flow  C                                    is "pipelined" in the vertical
150  C                              transport  C                                    so we need an fVer for each
151  C                              o uTrans: Zonal transport  C                                    variable.
152  C                              o vTrans: Meridional transport  C     phiHydC    :: hydrostatic potential anomaly at cell center
153  C                              o rTrans: Vertical transport  C                   In z coords phiHyd is the hydrostatic potential
154  C     maskC,maskUp             o maskC: land/water mask for tracer cells  C                      (=pressure/rho0) anomaly
155  C                              o maskUp: land/water mask for W points  C                   In p coords phiHyd is the geopotential height anomaly.
156  C     fVer[STUV]               o fVer: Vertical flux term - note fVer  C     phiHydF    :: hydrostatic potential anomaly at middle between 2 centers
157  C                                      is "pipelined" in the vertical  C     dPhiHydX,Y :: Gradient (X & Y directions) of hydrostatic potential anom.
158  C                                      so we need an fVer for each  C     phiSurfX,  ::  gradient of Surface potential (Pressure/rho, ocean)
159  C                                      variable.  C     phiSurfY             or geopotential (atmos) in X and Y direction
160  C     rhoK, rhoKM1   - Density at current level, and level above  C     guDissip   :: dissipation tendency (all explicit terms), u component
161  C     phiHyd         - Hydrostatic part of the potential phiHydi.  C     gvDissip   :: dissipation tendency (all explicit terms), v component
 C                      In z coords phiHydiHyd is the hydrostatic  
 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).  
162  C     iMin, iMax     - Ranges and sub-block indices on which calculations  C     iMin, iMax     - Ranges and sub-block indices on which calculations
163  C     jMin, jMax       are applied.  C     jMin, jMax       are applied.
164  C     bi, bj  C     bi, bj
165  C     k, kup,        - Index for layer above and below. kup and kDown  C     k, kup,        - Index for layer above and below. kup and kDown
166  C     kDown, km1       are switched with layer to be the appropriate  C     kDown, km1       are switched with layer to be the appropriate
167  C                      index into fVerTerm.  C                      index into fVerTerm.
       _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 maskC   (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)  
168        _RL fVerU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
169        _RL fVerV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
170        _RL phiHyd  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL phiHydF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
171        _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiHydC (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
172        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL dPhiHydX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
173          _RL dPhiHydY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
174        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
175        _RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
176        _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL guDissip(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
177        _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL gvDissip(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
178        _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
179        _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)  
   
 C This is currently used by IVDC and Diagnostics  
       _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
180    
181        INTEGER iMin, iMax        INTEGER iMin, iMax
182        INTEGER jMin, jMax        INTEGER jMin, jMax
183        INTEGER bi, bj        INTEGER bi, bj
184        INTEGER i, j        INTEGER i, j
185        INTEGER k, km1, kup, kDown        INTEGER k, km1, kp1, kup, kDown
186    
187    #ifdef ALLOW_DIAGNOSTICS
188          _RL tmpFac
189    #endif /* ALLOW_DIAGNOSTICS */
190    
 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)  
191    
192  C---    The algorithm...  C---    The algorithm...
193  C  C
# Line 162  C         salt* = salt[n] + dt x ( 3/2 G Line 232  C         salt* = salt[n] + dt x ( 3/2 G
232  C         (1 + dt * K * d_zz) theta[n] = theta*  C         (1 + dt * K * d_zz) theta[n] = theta*
233  C         (1 + dt * K * d_zz) salt[n] = salt*  C         (1 + dt * K * d_zz) salt[n] = salt*
234  C---  C---
235    CEOP
236    
237  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_DEBUG
238  C--   dummy statement to end declaration part        IF ( debugLevel .GE. debLevB )
239        ikey = 1       &   CALL DEBUG_ENTER( 'DYNAMICS', myThid )
240  #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  
         maskC  (i,j) = 0. _d 0  
         phiSurfX(i,j) = 0. _d 0  
         phiSurfY(i,j) = 0. _d 0  
        ENDDO  
       ENDDO  
241    
242    C-- Call to routine for calculation of
243    C   Eliassen-Palm-flux-forced U-tendency,
244    C   if desired:
245    #ifdef INCLUDE_EP_FORCING_CODE
246          CALL CALC_EP_FORCING(myThid)
247    #endif
248    
249  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
250  C--   HPF directive to help TAMC  C--   HPF directive to help TAMC
# Line 205  CHPF$ INDEPENDENT Line 255  CHPF$ INDEPENDENT
255    
256  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
257  C--    HPF directive to help TAMC  C--    HPF directive to help TAMC
258  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS,fVerU,fVerV  CHPF$  INDEPENDENT, NEW (fVerU,fVerV
259  CHPF$&                  ,phiHyd,utrans,vtrans,maskc,xA,yA  CHPF$&                  ,phiHydF
260  CHPF$&                  ,KappaRT,KappaRS,KappaRU,KappaRV  CHPF$&                  ,KappaRU,KappaRV
261  CHPF$&                  )  CHPF$&                  )
262  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
263    
# Line 216  CHPF$&                  ) Line 266  CHPF$&                  )
266  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
267            act1 = bi - myBxLo(myThid)            act1 = bi - myBxLo(myThid)
268            max1 = myBxHi(myThid) - myBxLo(myThid) + 1            max1 = myBxHi(myThid) - myBxLo(myThid) + 1
   
269            act2 = bj - myByLo(myThid)            act2 = bj - myByLo(myThid)
270            max2 = myByHi(myThid) - myByLo(myThid) + 1            max2 = myByHi(myThid) - myByLo(myThid) + 1
   
271            act3 = myThid - 1            act3 = myThid - 1
272            max3 = nTx*nTy            max3 = nTx*nTy
   
273            act4 = ikey_dynamics - 1            act4 = ikey_dynamics - 1
274              idynkey = (act1 + 1) + act2*max1
           ikey = (act1 + 1) + act2*max1  
275       &                      + act3*max1*max2       &                      + act3*max1*max2
276       &                      + act4*max1*max2*max3       &                      + act4*max1*max2*max3
277  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
278    
279  C--     Set up work arrays that need valid initial values  C--   Set up work arrays with valid (i.e. not NaN) values
280          DO j=1-OLy,sNy+OLy  C     These inital values do not alter the numerical results. They
281           DO i=1-OLx,sNx+OLx  C     just ensure that all memory references are to valid floating
282            rTrans(i,j)   = 0. _d 0  C     point numbers. This prevents spurious hardware signals due to
283            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  
284    
285          DO k=1,Nr          DO k=1,Nr
286           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
287            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
288  C This is currently also used by IVDC and Diagnostics             KappaRU(i,j,k) = 0. _d 0
289             ConvectCount(i,j,k) = 0.             KappaRV(i,j,k) = 0. _d 0
290             KappaRT(i,j,k) = 0. _d 0  #ifdef ALLOW_AUTODIFF_TAMC
291             KappaRS(i,j,k) = 0. _d 0  cph(
292    c--   need some re-initialisation here to break dependencies
293    cph)
294               gU(i,j,k,bi,bj) = 0. _d 0
295               gV(i,j,k,bi,bj) = 0. _d 0
296    #endif
297            ENDDO            ENDDO
298           ENDDO           ENDDO
299          ENDDO          ENDDO
300            DO j=1-OLy,sNy+OLy
301          iMin = 1-OLx+1           DO i=1-OLx,sNx+OLx
302          iMax = sNx+OLx            fVerU  (i,j,1) = 0. _d 0
303          jMin = 1-OLy+1            fVerU  (i,j,2) = 0. _d 0
304          jMax = sNy+OLy            fVerV  (i,j,1) = 0. _d 0
305              fVerV  (i,j,2) = 0. _d 0
306              phiHydF (i,j)  = 0. _d 0
307  #ifdef ALLOW_AUTODIFF_TAMC            phiHydC (i,j)  = 0. _d 0
308  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte            dPhiHydX(i,j)  = 0. _d 0
309  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte            dPhiHydY(i,j)  = 0. _d 0
310  CADJ STORE uvel(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte            phiSurfX(i,j)  = 0. _d 0
311  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte            phiSurfY(i,j)  = 0. _d 0
312  #endif /* ALLOW_AUTODIFF_TAMC */            guDissip(i,j)  = 0. _d 0
313              gvDissip(i,j)  = 0. _d 0
314  C--     Start of diagnostic loop           ENDDO
         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  
             CALL FIND_RHO(  
      I        bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,  
      I        theta, salt,  
      O        rhoK,  
      I        myThid )  
             IF (k.GT.1) CALL FIND_RHO(  
      I        bi, bj, iMin, iMax, jMin, jMax, k-1, k, eosType,  
      I        theta, salt,  
      O        rhoKm1,  
      I        myThid )  
             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)  
315          ENDDO          ENDDO
316    
317  #ifdef  ALLOW_OBCS  C--     Start computation of dynamics
318  C--     Calculate future values on open boundaries          iMin = 0
319          IF (useOBCS) THEN          iMax = sNx+1
320            CALL OBCS_CALC( bi, bj, myTime+deltaT,          jMin = 0
321       I            uVel, vVel, wVel, theta, salt,          jMax = sNy+1
      I            myThid )  
         ENDIF  
 #endif  /* ALLOW_OBCS */  
322    
 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_GMREDI  
 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  
323  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
324          ELSE  CADJ STORE wvel (:,:,:,bi,bj) =
325            DO k=1, Nr  CADJ &     comlev1_bibj, key = idynkey, byte = isbyte
             CALL GMREDI_CALC_TENSOR_DUMMY(  
      I             bi, bj, iMin, iMax, jMin, jMax, k,  
      I             sigmaX, sigmaY, sigmaR,  
      I             myThid )  
           ENDDO  
326  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
         ENDIF  
 #endif  /* ALLOW_GMREDI */  
327    
328  #ifdef  ALLOW_KPP  C--     Explicit part of the Surface Potentiel Gradient (add in TIMESTEP)
329  C--     Compute KPP mixing coefficients  C       (note: this loop will be replaced by CALL CALC_GRAD_ETA)
330          IF (useKPP) THEN          IF (implicSurfPress.NE.1.) THEN
331            CALL KPP_CALC(            CALL CALC_GRAD_PHI_SURF(
332       I                  bi, bj, myTime, myThid )       I         bi,bj,iMin,iMax,jMin,jMax,
333  #ifdef ALLOW_AUTODIFF_TAMC       I         etaN,
334          ELSE       O         phiSurfX,phiSurfY,
335            DO j=1-OLy,sNy+OLy       I         myThid )                        
             DO i=1-OLx,sNx+OLx  
               KPPhbl (i,j,bi,bj) = 1.0  
               KPPfrac(i,j,bi,bj) = 0.0  
               DO k = 1,Nr  
                  KPPghat   (i,j,k,bi,bj) = 0.0  
                  KPPviscAz (i,j,k,bi,bj) = viscAz  
                  KPPdiffKzT(i,j,k,bi,bj) = diffKzT  
                  KPPdiffKzS(i,j,k,bi,bj) = diffKzS  
               ENDDO  
             ENDDO  
           ENDDO  
 #endif /* ALLOW_AUTODIFF_TAMC */  
336          ENDIF          ENDIF
337    
338  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
339  CADJ STORE KPPghat   (:,:,:,bi,bj)  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=idynkey, byte=isbyte
340  CADJ &   , KPPviscAz (:,:,:,bi,bj)  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=idynkey, byte=isbyte
341  CADJ &   , KPPdiffKzT(:,:,:,bi,bj)  #ifdef ALLOW_KPP
342  CADJ &   , KPPdiffKzS(:,:,:,bi,bj)  CADJ STORE KPPviscAz (:,:,:,bi,bj)
343  CADJ &   , KPPfrac   (:,:  ,bi,bj)  CADJ &                 = comlev1_bibj, key=idynkey, byte=isbyte
344  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  
   
 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  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 C? Patrick Is this formula correct?  
 cph Yes, but I rewrote it.  
 cph Also, the KappaR? need the index k!  
          kkey = (ikey-1)*Nr + k  
 CADJ STORE KappaRT(:,:,k)    = comlev1_bibj_k, key = kkey, byte = isbyte  
 CADJ STORE KappaRS(:,:,k)    = comlev1_bibj_k, key = kkey, byte = isbyte  
345  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
346    
 C--      Get temporary terms used by tendency routines  
          CALL CALC_COMMON_FACTORS (  
      I        bi,bj,iMin,iMax,jMin,jMax,k,km1,kup,kDown,  
      O        xA,yA,uTrans,vTrans,rTrans,maskC,maskUp,  
      I        myThid)  
   
347  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL
348  C--      Calculate the total vertical diffusivity  C--      Calculate the total vertical diffusivity
349           CALL CALC_DIFFUSIVITY(          DO k=1,Nr
350             CALL CALC_VISCOSITY(
351       I        bi,bj,iMin,iMax,jMin,jMax,k,       I        bi,bj,iMin,iMax,jMin,jMax,k,
352       I        maskC,maskup,       O        KappaRU,KappaRV,
      O        KappaRT,KappaRS,KappaRU,KappaRV,  
353       I        myThid)       I        myThid)
354           ENDDO
355  #endif  #endif
356    
 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,maskC,  
      I         KappaRT,  
      U         fVerT,  
      I         myTime, myThid)  
            CALL TIMESTEP_TRACER(  
      I         bi,bj,iMin,iMax,jMin,jMax,k,  
      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,maskC,  
      I         KappaRS,  
      U         fVerS,  
      I         myTime, myThid)  
            CALL TIMESTEP_TRACER(  
      I         bi,bj,iMin,iMax,jMin,jMax,k,  
      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  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE gTNm1(:,:,k,bi,bj) = comlev1_bibj_k  
 CADJ &   , key = kkey, byte = isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
             CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )  
          END IF  
   
 C--     end of thermodynamic k loop (Nr:1)  
         ENDDO  
   
   
 #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  
357  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
358           idkey = iikey + 2  CADJ STORE KappaRU(:,:,:)
359  CADJ STORE gSNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ &                 = comlev1_bibj, key=idynkey, byte=isbyte
360    CADJ STORE KappaRV(:,:,:)
361    CADJ &                 = comlev1_bibj, key=idynkey, byte=isbyte
362  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
             CALL IMPLDIFF(  
      I         bi, bj, iMin, iMax, jMin, jMax,  
      I         deltaTtracer, KappaRS, recip_HFacC,  
      U         gSNm1,  
      I         myThid )  
          ENDIF  
   
 #ifdef   ALLOW_OBCS  
 C--      Apply open boundary conditions  
          IF (useOBCS) THEN  
            DO K=1,Nr  
              CALL OBCS_APPLY_TS( bi, bj, k, gTnm1, gSnm1, myThid )  
            ENDDO  
          END IF  
 #endif   /* ALLOW_OBCS */  
   
 C--     End If implicitDiffusion  
         ENDIF  
   
 C--     Start computation of dynamics  
         iMin = 1-OLx+2  
         iMax = sNx+OLx-1  
         jMin = 1-OLy+2  
         jMax = sNy+OLy-1  
   
 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  
363    
364  C--     Start of dynamics loop  C--     Start of dynamics loop
365          DO k=1,Nr          DO k=1,Nr
# Line 582  C--       kup    Cycles through 1,2 to p Line 369  C--       kup    Cycles through 1,2 to p
369  C--       kDown  Cycles through 2,1 to point to current layer  C--       kDown  Cycles through 2,1 to point to current layer
370    
371            km1  = MAX(1,k-1)            km1  = MAX(1,k-1)
372              kp1  = MIN(k+1,Nr)
373            kup  = 1+MOD(k+1,2)            kup  = 1+MOD(k+1,2)
374            kDown= 1+MOD(k,2)            kDown= 1+MOD(k,2)
375    
376    #ifdef ALLOW_AUTODIFF_TAMC
377             kkey = (idynkey-1)*Nr + k
378    c
379    CADJ STORE totphihyd (:,:,k,bi,bj)
380    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
381    CADJ STORE theta (:,:,k,bi,bj)
382    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
383    CADJ STORE salt  (:,:,k,bi,bj)
384    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
385    CADJ STORE gt(:,:,k,bi,bj)
386    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
387    CADJ STORE gs(:,:,k,bi,bj)
388    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
389    # ifdef NONLIN_FRSURF
390    cph-test
391    CADJ STORE  phiHydC (:,:)
392    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
393    CADJ STORE  phiHydF (:,:)
394    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
395    CADJ STORE  gudissip (:,:)
396    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
397    CADJ STORE  gvdissip (:,:)
398    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
399    CADJ STORE  fVerU (:,:,:)
400    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
401    CADJ STORE  fVerV (:,:,:)
402    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
403    CADJ STORE gu(:,:,k,bi,bj)
404    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
405    CADJ STORE gv(:,:,k,bi,bj)
406    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
407    CADJ STORE gunm1(:,:,k,bi,bj)
408    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
409    CADJ STORE gvnm1(:,:,k,bi,bj)
410    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
411    #  ifdef ALLOW_CD_CODE
412    CADJ STORE unm1(:,:,k,bi,bj)
413    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
414    CADJ STORE vnm1(:,:,k,bi,bj)
415    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
416    CADJ STORE uVelD(:,:,k,bi,bj)
417    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
418    CADJ STORE vVelD(:,:,k,bi,bj)
419    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
420    #  endif
421    # endif
422    #endif /* ALLOW_AUTODIFF_TAMC */
423    
424  C--      Integrate hydrostatic balance for phiHyd with BC of  C--      Integrate hydrostatic balance for phiHyd with BC of
425  C        phiHyd(z=0)=0  C        phiHyd(z=0)=0
426  C        distinguishe between Stagger and Non Stagger time stepping           IF ( implicitIntGravWave ) THEN
          IF (staggerTimeStep) THEN  
427             CALL CALC_PHI_HYD(             CALL CALC_PHI_HYD(
428       I        bi,bj,iMin,iMax,jMin,jMax,k,       I        bi,bj,iMin,iMax,jMin,jMax,k,
429       I        gTnm1, gSnm1,       I        gT, gS,
430       U        phiHyd,       U        phiHydF,
431       I        myThid )       O        phiHydC, dPhiHydX, dPhiHydY,
432         I        myTime, myIter, myThid )
433           ELSE           ELSE
434             CALL CALC_PHI_HYD(             CALL CALC_PHI_HYD(
435       I        bi,bj,iMin,iMax,jMin,jMax,k,       I        bi,bj,iMin,iMax,jMin,jMax,k,
436       I        theta, salt,       I        theta, salt,
437       U        phiHyd,       U        phiHydF,
438       I        myThid )       O        phiHydC, dPhiHydX, dPhiHydY,
439         I        myTime, myIter, myThid )
440           ENDIF           ENDIF
441    
442  C--      Calculate accelerations in the momentum equations (gU, gV, ...)  C--      Calculate accelerations in the momentum equations (gU, gV, ...)
443  C        and step forward storing the result in gUnm1, gVnm1, etc...  C        and step forward storing the result in gU, gV, etc...
444           IF ( momStepping ) THEN           IF ( momStepping ) THEN
445             CALL CALC_MOM_RHS(  #ifdef ALLOW_MOM_FLUXFORM
446               IF (.NOT. vectorInvariantMomentum) CALL MOM_FLUXFORM(
447         I         bi,bj,iMin,iMax,jMin,jMax,k,kup,kDown,
448         I         KappaRU, KappaRV,
449         U         fVerU, fVerV,
450         O         guDissip, gvDissip,
451         I         myTime, myIter, myThid)
452    #endif
453    #ifdef ALLOW_MOM_VECINV
454               IF (vectorInvariantMomentum) THEN
455    C
456    # ifdef ALLOW_AUTODIFF_TAMC
457    #  ifdef NONLIN_FRSURF
458    CADJ STORE fVerU(:,:,:)
459    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
460    CADJ STORE fVerV(:,:,:)
461    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
462    #  endif
463    # endif /* ALLOW_AUTODIFF_TAMC */
464    C
465                 CALL MOM_VECINV(
466       I         bi,bj,iMin,iMax,jMin,jMax,k,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax,k,kup,kDown,
467       I         phiHyd,KappaRU,KappaRV,       I         KappaRU, KappaRV,
468       U         fVerU, fVerV,       U         fVerU, fVerV,
469       I         myTime, myThid)       O         guDissip, gvDissip,
470         I         myTime, myIter, myThid)
471               ENDIF
472    #endif
473             CALL TIMESTEP(             CALL TIMESTEP(
474       I         bi,bj,iMin,iMax,jMin,jMax,k,       I         bi,bj,iMin,iMax,jMin,jMax,k,
475       I         phiHyd, phiSurfX, phiSurfY,       I         dPhiHydX,dPhiHydY, phiSurfX, phiSurfY,
476       I         myIter, myThid)       I         guDissip, gvDissip,
477         I         myTime, myIter, myThid)
478    
479  #ifdef   ALLOW_OBCS  #ifdef   ALLOW_OBCS
480  C--      Apply open boundary conditions  C--      Apply open boundary conditions
481           IF (useOBCS) THEN             IF (useOBCS) THEN
482             CALL OBCS_APPLY_UV( bi, bj, k, gUnm1, gVnm1, myThid )               CALL OBCS_APPLY_UV( bi, bj, k, gU, gV, myThid )
483           END IF             ENDIF
484  #endif   /* ALLOW_OBCS */  #endif   /* ALLOW_OBCS */
485    
 #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 */  
486           ENDIF           ENDIF
487    
488    
489  C--     end of dynamics k loop (1:Nr)  C--     end of dynamics k loop (1:Nr)
490          ENDDO          ENDDO
491    
492    C--     Implicit Vertical advection & viscosity
493    #if (defined (INCLUDE_IMPLVERTADV_CODE) && defined (ALLOW_MOM_COMMON))
494  C--     Implicit viscosity          IF ( momImplVertAdv ) THEN
495          IF (implicitViscosity.AND.momStepping) THEN            CALL MOM_U_IMPLICIT_R( kappaRU,
496         I                           bi, bj, myTime, myIter, myThid )
497              CALL MOM_V_IMPLICIT_R( kappaRV,
498         I                           bi, bj, myTime, myIter, myThid )
499            ELSEIF ( implicitViscosity ) THEN
500    #else /* INCLUDE_IMPLVERTADV_CODE */
501            IF     ( implicitViscosity ) THEN
502    #endif /* INCLUDE_IMPLVERTADV_CODE */
503  #ifdef    ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
504            idkey = iikey + 3  CADJ STORE KappaRU(:,:,:) = comlev1_bibj , key=idynkey, byte=isbyte
505  CADJ STORE gUNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE gU(:,:,:,bi,bj) = comlev1_bibj , key=idynkey, byte=isbyte
506  #endif    /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
507            CALL IMPLDIFF(            CALL IMPLDIFF(
508       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
509       I         deltaTmom, KappaRU,recip_HFacW,       I         -1, KappaRU,recip_HFacW,
510       U         gUNm1,       U         gU,
511       I         myThid )       I         myThid )
512  #ifdef    ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
513            idkey = iikey + 4  CADJ STORE KappaRV(:,:,:) = comlev1_bibj , key=idynkey, byte=isbyte
514  CADJ STORE gVNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE gV(:,:,:,bi,bj) = comlev1_bibj , key=idynkey, byte=isbyte
515  #endif    /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
516            CALL IMPLDIFF(            CALL IMPLDIFF(
517       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
518       I         deltaTmom, KappaRV,recip_HFacS,       I         -2, KappaRV,recip_HFacS,
519       U         gVNm1,       U         gV,
520       I         myThid )       I         myThid )
521            ENDIF
522    
523  #ifdef   ALLOW_OBCS  #ifdef   ALLOW_OBCS
524  C--      Apply open boundary conditions  C--      Apply open boundary conditions
525           IF (useOBCS) THEN          IF ( useOBCS .AND.(implicitViscosity.OR.momImplVertAdv) ) THEN
526             DO K=1,Nr             DO K=1,Nr
527               CALL OBCS_APPLY_UV( bi, bj, k, gUnm1, gVnm1, myThid )               CALL OBCS_APPLY_UV( bi, bj, k, gU, gV, myThid )
528             ENDDO             ENDDO
529           END IF          ENDIF
530  #endif   /* ALLOW_OBCS */  #endif   /* ALLOW_OBCS */
531    
532  #ifdef    INCLUDE_CD_CODE  #ifdef    ALLOW_CD_CODE
533            IF (implicitViscosity.AND.useCDscheme) THEN
534  #ifdef    ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
535            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  
536  #endif    /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
537            CALL IMPLDIFF(            CALL IMPLDIFF(
538       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
539       I         deltaTmom, KappaRU,recip_HFacW,       I         0, KappaRU,recip_HFacW,
540       U         vVelD,       U         vVelD,
541       I         myThid )       I         myThid )
542  #ifdef    ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
543            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  
544  #endif    /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
545            CALL IMPLDIFF(            CALL IMPLDIFF(
546       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
547       I         deltaTmom, KappaRV,recip_HFacS,       I         0, KappaRV,recip_HFacS,
548       U         uVelD,       U         uVelD,
549       I         myThid )       I         myThid )
 #endif    /* INCLUDE_CD_CODE */  
 C--     End If implicitViscosity.AND.momStepping  
550          ENDIF          ENDIF
551    #endif    /* ALLOW_CD_CODE */
552    C--     End implicit Vertical advection & viscosity
553    
 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_CUMULATE(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 */  
   
554         ENDDO         ENDDO
555        ENDDO        ENDDO
556    
557    #ifdef ALLOW_OBCS
558          IF (useOBCS) THEN
559           CALL OBCS_PRESCRIBE_EXCHANGES(myThid)
560          ENDIF
561    #endif
562    
563    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
564    
565    #ifdef ALLOW_NONHYDROSTATIC
566    C--   Step forward W field in N-H algorithm
567          IF ( nonHydrostatic ) THEN
568    #ifdef ALLOW_DEBUG
569             IF ( debugLevel .GE. debLevB )
570         &     CALL DEBUG_CALL('CALC_GW', myThid )
571    #endif
572             CALL TIMER_START('CALC_GW          [DYNAMICS]',myThid)
573             CALL CALC_GW( myTime, myIter, myThid )
574          ENDIF
575          IF ( nonHydrostatic.OR.implicitIntGravWave )
576         &   CALL TIMESTEP_WVEL( myTime, myIter, myThid )
577          IF ( nonHydrostatic )
578         &   CALL TIMER_STOP ('CALC_GW          [DYNAMICS]',myThid)
579    #endif
580    
581    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
582    
583    Cml(
584    C     In order to compare the variance of phiHydLow of a p/z-coordinate
585    C     run with etaH of a z/p-coordinate run the drift of phiHydLow
586    C     has to be removed by something like the following subroutine:
587    C      CALL REMOVE_MEAN_RL( 1, phiHydLow, maskH, maskH, rA, drF,
588    C     &                'phiHydLow', myThid )
589    Cml)
590    
591    #ifdef ALLOW_DIAGNOSTICS
592          IF ( useDiagnostics ) THEN
593    
594           CALL DIAGNOSTICS_FILL(totPhihyd,'PHIHYD  ',0,Nr,0,1,1,myThid)
595           CALL DIAGNOSTICS_FILL(phiHydLow,'PHIBOT  ',0, 1,0,1,1,myThid)
596    
597           tmpFac = 1. _d 0
598           CALL DIAGNOSTICS_SCALE_FILL(totPhihyd,tmpFac,2,
599         &                                 'PHIHYDSQ',0,Nr,0,1,1,myThid)
600    
601           CALL DIAGNOSTICS_SCALE_FILL(phiHydLow,tmpFac,2,
602         &                                 'PHIBOTSQ',0, 1,0,1,1,myThid)
603    
604          ENDIF
605    #endif /* ALLOW_DIAGNOSTICS */
606          
607    #ifdef ALLOW_DEBUG
608          If ( debugLevel .GE. debLevB ) THEN
609           CALL DEBUG_STATS_RL(1,EtaN,'EtaN (DYNAMICS)',myThid)
610           CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (DYNAMICS)',myThid)
611           CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (DYNAMICS)',myThid)
612           CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (DYNAMICS)',myThid)
613           CALL DEBUG_STATS_RL(Nr,theta,'Theta (DYNAMICS)',myThid)
614           CALL DEBUG_STATS_RL(Nr,salt,'Salt (DYNAMICS)',myThid)
615           CALL DEBUG_STATS_RL(Nr,gU,'Gu (DYNAMICS)',myThid)
616           CALL DEBUG_STATS_RL(Nr,gV,'Gv (DYNAMICS)',myThid)
617           CALL DEBUG_STATS_RL(Nr,gT,'Gt (DYNAMICS)',myThid)
618           CALL DEBUG_STATS_RL(Nr,gS,'Gs (DYNAMICS)',myThid)
619    #ifndef ALLOW_ADAMSBASHFORTH_3
620           CALL DEBUG_STATS_RL(Nr,guNm1,'GuNm1 (DYNAMICS)',myThid)
621           CALL DEBUG_STATS_RL(Nr,gvNm1,'GvNm1 (DYNAMICS)',myThid)
622           CALL DEBUG_STATS_RL(Nr,gtNm1,'GtNm1 (DYNAMICS)',myThid)
623           CALL DEBUG_STATS_RL(Nr,gsNm1,'GsNm1 (DYNAMICS)',myThid)
624    #endif
625          ENDIF
626    #endif
627    
628    #ifdef DYNAMICS_GUGV_EXCH_CHECK
629    C- jmc: For safety checking only: This Exchange here should not change
630    C       the solution. If solution changes, it means something is wrong,
631    C       but it does not mean that it is less wrong with this exchange.
632          IF ( debugLevel .GT. debLevB ) THEN
633           CALL EXCH_UV_XYZ_RL(gU,gV,.TRUE.,myThid)
634          ENDIF
635    #endif
636    
637    #ifdef ALLOW_DEBUG
638          IF ( debugLevel .GE. debLevB )
639         &   CALL DEBUG_LEAVE( 'DYNAMICS', myThid )
640    #endif
641    
642        RETURN        RETURN
643        END        END
Legend:
Removed from v.1.66  
changed lines
  Added in v.1.131
  ViewVC Help
Powered by ViewVC 1.1.22