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