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

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

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

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