/[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.146 by jmc, Fri May 14 23:21:02 2010 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     !FUNCTIONS:
150    #ifdef ALLOW_DIAGNOSTICS
151          LOGICAL  DIAGNOSTICS_IS_ON
152          EXTERNAL DIAGNOSTICS_IS_ON
153    #endif
154    
155    C     !LOCAL VARIABLES:
156  C     == Local variables  C     == Local variables
157  C     xA, yA                 - Per block temporaries holding face areas  C     fVer[UV]               o fVer: Vertical flux term - note fVer
158  C     uTrans, vTrans, rTrans - Per block temporaries holding flow  C                                    is "pipelined" in the vertical
159  C                              transport  C                                    so we need an fVer for each
160  C                              o uTrans: Zonal transport  C                                    variable.
161  C                              o vTrans: Meridional transport  C     phiHydC    :: hydrostatic potential anomaly at cell center
162  C                              o rTrans: Vertical transport  C                   In z coords phiHyd is the hydrostatic potential
163  C     maskUp                   o maskUp: land/water mask for W points  C                      (=pressure/rho0) anomaly
164  C     fVer[STUV]               o fVer: Vertical flux term - note fVer  C                   In p coords phiHyd is the geopotential height anomaly.
165  C                                      is "pipelined" in the vertical  C     phiHydF    :: hydrostatic potential anomaly at middle between 2 centers
166  C                                      so we need an fVer for each  C     dPhiHydX,Y :: Gradient (X & Y directions) of hydrostatic potential anom.
167  C                                      variable.  C     phiSurfX,  ::  gradient of Surface potential (Pressure/rho, ocean)
168  C     rhoK, rhoKM1   - Density at current level, and level above  C     phiSurfY             or geopotential (atmos) in X and Y direction
169  C     phiHyd         - Hydrostatic part of the potential phiHydi.  C     guDissip   :: dissipation tendency (all explicit terms), u component
170  C                      In z coords phiHydiHyd is the hydrostatic  C     gvDissip   :: dissipation tendency (all explicit terms), v component
171  C                      Potential (=pressure/rho0) anomaly  C     KappaRU    :: vertical viscosity
172  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).  
173  C     iMin, iMax     - Ranges and sub-block indices on which calculations  C     iMin, iMax     - Ranges and sub-block indices on which calculations
174  C     jMin, jMax       are applied.  C     jMin, jMax       are applied.
175  C     bi, bj  C     bi, bj
176  C     k, kup,        - Index for layer above and below. kup and kDown  C     k, kup,        - Index for layer above and below. kup and kDown
177  C     kDown, km1       are switched with layer to be the appropriate  C     kDown, km1       are switched with layer to be the appropriate
178  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)  
179        _RL fVerU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
180        _RL fVerV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
181        _RL phiHyd  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL phiHydF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
182        _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiHydC (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
183        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL dPhiHydX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
184          _RL dPhiHydY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
185        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
186        _RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
187        _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL guDissip(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
188        _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL gvDissip(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
189        _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
190        _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)  
191    
192        INTEGER iMin, iMax        INTEGER iMin, iMax
193        INTEGER jMin, jMax        INTEGER jMin, jMax
194        INTEGER bi, bj        INTEGER bi, bj
195        INTEGER i, j        INTEGER i, j
196        INTEGER k, km1, kup, kDown        INTEGER k, km1, kp1, kup, kDown
197    
198    #ifdef ALLOW_DIAGNOSTICS
199          LOGICAL dPhiHydDiagIsOn
200          _RL tmpFac
201    #endif /* ALLOW_DIAGNOSTICS */
202    
203    
 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)  
   
204  C---    The algorithm...  C---    The algorithm...
205  C  C
206  C       "Correction Step"  C       "Correction Step"
# Line 165  C         salt* = salt[n] + dt x ( 3/2 G Line 244  C         salt* = salt[n] + dt x ( 3/2 G
244  C         (1 + dt * K * d_zz) theta[n] = theta*  C         (1 + dt * K * d_zz) theta[n] = theta*
245  C         (1 + dt * K * d_zz) salt[n] = salt*  C         (1 + dt * K * d_zz) salt[n] = salt*
246  C---  C---
247    CEOP
248    
249  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_DEBUG
250  C--   dummy statement to end declaration part        IF ( debugLevel .GE. debLevB )
251        ikey = 1       &   CALL DEBUG_ENTER( 'DYNAMICS', myThid )
252  #endif /* ALLOW_AUTODIFF_TAMC */  #endif
253    
254  C--   Set up work arrays with valid (i.e. not NaN) values  #ifdef ALLOW_DIAGNOSTICS
255  C     These inital values do not alter the numerical results. They        dPhiHydDiagIsOn = .FALSE.
256  C     just ensure that all memory references are to valid floating        IF ( useDiagnostics )
257  C     point numbers. This prevents spurious hardware signals due to       &  dPhiHydDiagIsOn = DIAGNOSTICS_IS_ON( 'Um_dPHdx', myThid )
258  C     uninitialised but inert locations.       &               .OR. DIAGNOSTICS_IS_ON( 'Vm_dPHdy', myThid )
259        DO j=1-OLy,sNy+OLy  #endif
        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  
260    
261    C-- Call to routine for calculation of
262    C   Eliassen-Palm-flux-forced U-tendency,
263    C   if desired:
264    #ifdef INCLUDE_EP_FORCING_CODE
265          CALL CALC_EP_FORCING(myThid)
266    #endif
267    
268  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
269  C--   HPF directive to help TAMC  C--   HPF directive to help TAMC
# Line 207  CHPF$ INDEPENDENT Line 274  CHPF$ INDEPENDENT
274    
275  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
276  C--    HPF directive to help TAMC  C--    HPF directive to help TAMC
277  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS,fVerU,fVerV  CHPF$  INDEPENDENT, NEW (fVerU,fVerV
278  CHPF$&                  ,phiHyd,utrans,vtrans,xA,yA  CHPF$&                  ,phiHydF
279  CHPF$&                  ,KappaRT,KappaRS,KappaRU,KappaRV  CHPF$&                  ,KappaRU,KappaRV
280  CHPF$&                  )  CHPF$&                  )
281  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
282    
# Line 218  CHPF$&                  ) Line 285  CHPF$&                  )
285  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
286            act1 = bi - myBxLo(myThid)            act1 = bi - myBxLo(myThid)
287            max1 = myBxHi(myThid) - myBxLo(myThid) + 1            max1 = myBxHi(myThid) - myBxLo(myThid) + 1
   
288            act2 = bj - myByLo(myThid)            act2 = bj - myByLo(myThid)
289            max2 = myByHi(myThid) - myByLo(myThid) + 1            max2 = myByHi(myThid) - myByLo(myThid) + 1
   
290            act3 = myThid - 1            act3 = myThid - 1
291            max3 = nTx*nTy            max3 = nTx*nTy
   
292            act4 = ikey_dynamics - 1            act4 = ikey_dynamics - 1
293              idynkey = (act1 + 1) + act2*max1
           ikey = (act1 + 1) + act2*max1  
294       &                      + act3*max1*max2       &                      + act3*max1*max2
295       &                      + act4*max1*max2*max3       &                      + act4*max1*max2*max3
296  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
297    
298  C--     Set up work arrays that need valid initial values  C--   Set up work arrays with valid (i.e. not NaN) values
299          DO j=1-OLy,sNy+OLy  C     These inital values do not alter the numerical results. They
300           DO i=1-OLx,sNx+OLx  C     just ensure that all memory references are to valid floating
301            rTrans(i,j)   = 0. _d 0  C     point numbers. This prevents spurious hardware signals due to
302            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  
303    
304    #ifdef ALLOW_AUTODIFF_TAMC
305          DO k=1,Nr          DO k=1,Nr
306           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
307            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
308  C This is currently also used by IVDC and Diagnostics             KappaRU(i,j,k) = 0. _d 0
309             ConvectCount(i,j,k) = 0.             KappaRV(i,j,k) = 0. _d 0
310             KappaRT(i,j,k) = 0. _d 0  cph(
311             KappaRS(i,j,k) = 0. _d 0  c--   need some re-initialisation here to break dependencies
312    cph)
313               gU(i,j,k,bi,bj) = 0. _d 0
314               gV(i,j,k,bi,bj) = 0. _d 0
315            ENDDO            ENDDO
316           ENDDO           ENDDO
317          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  
318  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
319               CALL FIND_RHO(          DO j=1-OLy,sNy+OLy
320       I        bi, bj, iMin, iMax, jMin, jMax, k-1, k, eosType,           DO i=1-OLx,sNx+OLx
321       I        theta, salt,            fVerU  (i,j,1) = 0. _d 0
322       O        rhoKm1,            fVerU  (i,j,2) = 0. _d 0
323       I        myThid )            fVerV  (i,j,1) = 0. _d 0
324              ENDIF            fVerV  (i,j,2) = 0. _d 0
325              CALL GRAD_SIGMA(            phiHydF (i,j)  = 0. _d 0
326       I             bi, bj, iMin, iMax, jMin, jMax, k,            phiHydC (i,j)  = 0. _d 0
327       I             rhoK, rhoKm1, rhoK,  #ifndef INCLUDE_PHIHYD_CALCULATION_CODE
328       O             sigmaX, sigmaY, sigmaR,            dPhiHydX(i,j)  = 0. _d 0
329       I             myThid )            dPhiHydY(i,j)  = 0. _d 0
330            ENDIF  #endif
331              phiSurfX(i,j)  = 0. _d 0
332  C--       Implicit Vertical Diffusion for Convection            phiSurfY(i,j)  = 0. _d 0
333  c ==> should use sigmaR !!!            guDissip(i,j)  = 0. _d 0
334            IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN            gvDissip(i,j)  = 0. _d 0
335              CALL CALC_IVDC(  #ifdef ALLOW_AUTODIFF_TAMC
336       I        bi, bj, iMin, iMax, jMin, jMax, k,            phiHydLow(i,j,bi,bj) = 0. _d 0
337       I        rhoKm1, rhoK,  # ifdef NONLIN_FRSURF
338       U        ConvectCount, KappaRT, KappaRS,  #  ifndef DISABLE_RSTAR_CODE
339       I        myTime, myIter, myThid)            dWtransC(i,j,bi,bj) = 0. _d 0
340            ENDIF            dWtransU(i,j,bi,bj) = 0. _d 0
341              dWtransV(i,j,bi,bj) = 0. _d 0
342  C--     end of diagnostic k loop (Nr:1)  #  endif
343    # endif
344    #endif
345             ENDDO
346          ENDDO          ENDDO
347    
348  #ifdef ALLOW_AUTODIFF_TAMC  C--     Start computation of dynamics
349  cph avoids recomputation of integrate_for_w          iMin = 0
350  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte          iMax = sNx+1
351  #endif /* ALLOW_AUTODIFF_TAMC */          jMin = 0
352            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  
353    
354  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
355  CADJ STORE Kwx(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE wvel (:,:,:,bi,bj) =
356  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  
357  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
358    
359  #endif  /* ALLOW_GMREDI */  C--     Explicit part of the Surface Potentiel Gradient (add in TIMESTEP)
360    C       (note: this loop will be replaced by CALL CALC_GRAD_ETA)
361  #ifdef  ALLOW_KPP          IF (implicSurfPress.NE.1.) THEN
362  C--     Compute KPP mixing coefficients            CALL CALC_GRAD_PHI_SURF(
363          IF (useKPP) THEN       I         bi,bj,iMin,iMax,jMin,jMax,
364            CALL KPP_CALC(       I         etaN,
365       I                  bi, bj, myTime, myThid )       O         phiSurfX,phiSurfY,
366  #ifdef ALLOW_AUTODIFF_TAMC       I         myThid )
         ELSE  
           CALL KPP_CALC_DUMMY(  
      I                  bi, bj, myTime, myThid )  
 #endif /* ALLOW_AUTODIFF_TAMC */  
367          ENDIF          ENDIF
368    
369  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
370  CADJ STORE KPPghat   (:,:,:,bi,bj)  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=idynkey, byte=isbyte
371  CADJ &   , KPPviscAz (:,:,:,bi,bj)  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=idynkey, byte=isbyte
372  CADJ &   , KPPdiffKzT(:,:,:,bi,bj)  #ifdef ALLOW_KPP
373  CADJ &   , KPPdiffKzS(:,:,:,bi,bj)  CADJ STORE KPPviscAz (:,:,:,bi,bj)
374  CADJ &   , KPPfrac   (:,:  ,bi,bj)  CADJ &                 = comlev1_bibj, key=idynkey, byte=isbyte
375  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  
376  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
377    
378  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL
379  C--      Calculate the total vertical diffusivity  C--     Calculate the total vertical viscosity
380           CALL CALC_DIFFUSIVITY(          CALL CALC_VISCOSITY(
381       I        bi,bj,iMin,iMax,jMin,jMax,k,       I            bi,bj, iMin,iMax,jMin,jMax,
382       I        maskUp,       O            KappaRU, KappaRV,
383       O        KappaRT,KappaRS,KappaRU,KappaRV,       I            myThid )
384       I        myThid)  #else
385  #endif          DO k=1,Nr
386             DO j=1-OLy,sNy+OLy
387  C--      Calculate active tracer tendencies (gT,gS,...)            DO i=1-OLx,sNx+OLx
388  C        and step forward storing result in gTnm1, gSnm1, etc.             KappaRU(i,j,k) = 0. _d 0
389           IF ( tempStepping ) THEN             KappaRV(i,j,k) = 0. _d 0
390             CALL CALC_GT(            ENDDO
391       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)  
392          ENDDO          ENDDO
393    #endif
394    
395  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
396  C? Patrick? What about this one?  CADJ STORE KappaRU(:,:,:)
397  cph Keys iikey and idkey don't seem to be needed  CADJ &     = comlev1_bibj, key=idynkey, byte=isbyte
398  cph since storing occurs on different tape for each  CADJ STORE KappaRV(:,:,:)
399  cph impldiff call anyways.  CADJ &     = comlev1_bibj, key=idynkey, byte=isbyte
 cph Thus, common block comlev1_impl isn't needed either.  
 cph Storing below needed in the case useGMREDI.  
         iikey = (ikey-1)*maximpl  
400  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
401    
 C--     Implicit diffusion  
         IF (implicitDiffusion) THEN  
   
          IF (tempStepping) THEN  
 #ifdef ALLOW_AUTODIFF_TAMC  
             idkey = iikey + 1  
 CADJ STORE gTNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
             CALL IMPLDIFF(  
      I         bi, bj, iMin, iMax, jMin, jMax,  
      I         deltaTtracer, KappaRT, recip_HFacC,  
      U         gTNm1,  
      I         myThid )  
          ENDIF  
   
          IF (saltStepping) THEN  
 #ifdef ALLOW_AUTODIFF_TAMC  
          idkey = iikey + 2  
 CADJ STORE gSNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
             CALL IMPLDIFF(  
      I         bi, bj, iMin, iMax, jMin, jMax,  
      I         deltaTtracer, KappaRS, recip_HFacC,  
      U         gSNm1,  
      I         myThid )  
          ENDIF  
   
 #ifdef   ALLOW_OBCS  
 C--      Apply open boundary conditions  
          IF (useOBCS) THEN  
            DO K=1,Nr  
              CALL OBCS_APPLY_TS( bi, bj, k, gTnm1, gSnm1, myThid )  
            ENDDO  
          END IF  
 #endif   /* ALLOW_OBCS */  
   
 C--     End If implicitDiffusion  
         ENDIF  
   
 C--     Start computation of dynamics  
         iMin = 1-OLx+2  
         iMax = sNx+OLx-1  
         jMin = 1-OLy+2  
         jMax = sNy+OLy-1  
   
 C--     Explicit part of the Surface Potentiel Gradient (add in TIMESTEP)  
 C       (note: this loop will be replaced by CALL CALC_GRAD_ETA)  
         IF (implicSurfPress.NE.1.) THEN  
           CALL CALC_GRAD_PHI_SURF(  
      I         bi,bj,iMin,iMax,jMin,jMax,  
      I         etaN,  
      O         phiSurfX,phiSurfY,  
      I         myThid )                          
         ENDIF  
   
402  C--     Start of dynamics loop  C--     Start of dynamics loop
403          DO k=1,Nr          DO k=1,Nr
404    
# Line 617  C--       kup    Cycles through 1,2 to p Line 407  C--       kup    Cycles through 1,2 to p
407  C--       kDown  Cycles through 2,1 to point to current layer  C--       kDown  Cycles through 2,1 to point to current layer
408    
409            km1  = MAX(1,k-1)            km1  = MAX(1,k-1)
410              kp1  = MIN(k+1,Nr)
411            kup  = 1+MOD(k+1,2)            kup  = 1+MOD(k+1,2)
412            kDown= 1+MOD(k,2)            kDown= 1+MOD(k,2)
413    
414  C--      Integrate hydrostatic balance for phiHyd with BC of  #ifdef ALLOW_AUTODIFF_TAMC
415             kkey = (idynkey-1)*Nr + k
416    c
417    CADJ STORE totphihyd (:,:,k,bi,bj)
418    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
419    CADJ STORE phihydlow (:,:,bi,bj)
420    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
421    CADJ STORE theta (:,:,k,bi,bj)
422    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
423    CADJ STORE salt  (:,:,k,bi,bj)
424    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
425    CADJ STORE gt(:,:,k,bi,bj)
426    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
427    CADJ STORE gs(:,:,k,bi,bj)
428    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
429    # ifdef NONLIN_FRSURF
430    cph-test
431    CADJ STORE  phiHydC (:,:)
432    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
433    CADJ STORE  phiHydF (:,:)
434    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
435    CADJ STORE  gudissip (:,:)
436    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
437    CADJ STORE  gvdissip (:,:)
438    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
439    CADJ STORE  fVerU (:,:,:)
440    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
441    CADJ STORE  fVerV (:,:,:)
442    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
443    CADJ STORE gu(:,:,k,bi,bj)
444    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
445    CADJ STORE gv(:,:,k,bi,bj)
446    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
447    CADJ STORE gunm1(:,:,k,bi,bj)
448    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
449    CADJ STORE gvnm1(:,:,k,bi,bj)
450    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
451    #  ifdef ALLOW_CD_CODE
452    CADJ STORE unm1(:,:,k,bi,bj)
453    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
454    CADJ STORE vnm1(:,:,k,bi,bj)
455    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
456    CADJ STORE uVelD(:,:,k,bi,bj)
457    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
458    CADJ STORE vVelD(:,:,k,bi,bj)
459    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
460    #  endif
461    # endif
462    # ifdef ALLOW_DEPTH_CONTROL
463    CADJ STORE  fVerU (:,:,:)
464    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
465    CADJ STORE  fVerV (:,:,:)
466    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
467    # endif
468    #endif /* ALLOW_AUTODIFF_TAMC */
469    
470    C--      Integrate hydrostatic balance for phiHyd with BC of
471  C        phiHyd(z=0)=0  C        phiHyd(z=0)=0
472  C        distinguishe between Stagger and Non Stagger time stepping           IF ( implicitIntGravWave ) THEN
          IF (staggerTimeStep) THEN  
473             CALL CALC_PHI_HYD(             CALL CALC_PHI_HYD(
474       I        bi,bj,iMin,iMax,jMin,jMax,k,       I        bi,bj,iMin,iMax,jMin,jMax,k,
475       I        gTnm1, gSnm1,       I        gT, gS,
476       U        phiHyd,       U        phiHydF,
477       I        myThid )       O        phiHydC, dPhiHydX, dPhiHydY,
478         I        myTime, myIter, myThid )
479           ELSE           ELSE
480             CALL CALC_PHI_HYD(             CALL CALC_PHI_HYD(
481       I        bi,bj,iMin,iMax,jMin,jMax,k,       I        bi,bj,iMin,iMax,jMin,jMax,k,
482       I        theta, salt,       I        theta, salt,
483       U        phiHyd,       U        phiHydF,
484       I        myThid )       O        phiHydC, dPhiHydX, dPhiHydY,
485         I        myTime, myIter, myThid )
486             ENDIF
487    #ifdef ALLOW_DIAGNOSTICS
488             IF ( dPhiHydDiagIsOn ) THEN
489               tmpFac = -1. _d 0
490               CALL DIAGNOSTICS_SCALE_FILL( dPhiHydX, tmpFac, 1,
491         &                           'Um_dPHdx', k, 1, 2, bi, bj, myThid )
492               CALL DIAGNOSTICS_SCALE_FILL( dPhiHydY, tmpFac, 1,
493         &                           'Vm_dPHdy', k, 1, 2, bi, bj, myThid )
494           ENDIF           ENDIF
495    #endif /* ALLOW_DIAGNOSTICS */
496    
497  C--      Calculate accelerations in the momentum equations (gU, gV, ...)  C--      Calculate accelerations in the momentum equations (gU, gV, ...)
498  C        and step forward storing the result in gUnm1, gVnm1, etc...  C        and step forward storing the result in gU, gV, etc...
499           IF ( momStepping ) THEN           IF ( momStepping ) THEN
500             CALL CALC_MOM_RHS(  #ifdef ALLOW_AUTODIFF_TAMC
501    # ifdef NONLIN_FRSURF
502    #  ifndef DISABLE_RSTAR_CODE
503    CADJ STORE dWtransC(:,:,bi,bj)
504    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
505    CADJ STORE dWtransU(:,:,bi,bj)
506    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
507    CADJ STORE dWtransV(:,:,bi,bj)
508    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
509    #  endif
510    # endif
511    #endif
512               IF (.NOT. vectorInvariantMomentum) THEN
513    #ifdef ALLOW_MOM_FLUXFORM
514    C
515                  CALL MOM_FLUXFORM(
516         I         bi,bj,iMin,iMax,jMin,jMax,k,kup,kDown,
517         I         KappaRU, KappaRV,
518         U         fVerU, fVerV,
519         O         guDissip, gvDissip,
520         I         myTime, myIter, myThid)
521    #endif
522               ELSE
523    #ifdef ALLOW_MOM_VECINV
524    C
525    # ifdef ALLOW_AUTODIFF_TAMC
526    #  ifdef NONLIN_FRSURF
527    CADJ STORE fVerU(:,:,:)
528    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
529    CADJ STORE fVerV(:,:,:)
530    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
531    #  endif
532    # endif /* ALLOW_AUTODIFF_TAMC */
533    C
534                 CALL MOM_VECINV(
535       I         bi,bj,iMin,iMax,jMin,jMax,k,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax,k,kup,kDown,
536       I         phiHyd,KappaRU,KappaRV,       I         KappaRU, KappaRV,
537       U         fVerU, fVerV,       U         fVerU, fVerV,
538       I         myTime, myThid)       O         guDissip, gvDissip,
539         I         myTime, myIter, myThid)
540    #endif
541               ENDIF
542    C
543             CALL TIMESTEP(             CALL TIMESTEP(
544       I         bi,bj,iMin,iMax,jMin,jMax,k,       I         bi,bj,iMin,iMax,jMin,jMax,k,
545       I         phiHyd, phiSurfX, phiSurfY,       I         dPhiHydX,dPhiHydY, phiSurfX, phiSurfY,
546       I         myIter, myThid)       I         guDissip, gvDissip,
547         I         myTime, myIter, myThid)
548    
549  #ifdef   ALLOW_OBCS  #ifdef   ALLOW_OBCS
550  C--      Apply open boundary conditions  C--      Apply open boundary conditions
551           IF (useOBCS) THEN             IF (useOBCS) THEN
552             CALL OBCS_APPLY_UV( bi, bj, k, gUnm1, gVnm1, myThid )               CALL OBCS_APPLY_UV( bi, bj, k, gU, gV, myThid )
553           END IF             ENDIF
554  #endif   /* ALLOW_OBCS */  #endif   /* ALLOW_OBCS */
555    
 #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 */  
556           ENDIF           ENDIF
557    
558    
559  C--     end of dynamics k loop (1:Nr)  C--     end of dynamics k loop (1:Nr)
560          ENDDO          ENDDO
561    
562    C--     Implicit Vertical advection & viscosity
563    #if (defined (INCLUDE_IMPLVERTADV_CODE) && defined (ALLOW_MOM_COMMON))
564  C--     Implicit viscosity          IF ( momImplVertAdv ) THEN
565          IF (implicitViscosity.AND.momStepping) THEN            CALL MOM_U_IMPLICIT_R( kappaRU,
566         I                           bi, bj, myTime, myIter, myThid )
567              CALL MOM_V_IMPLICIT_R( kappaRV,
568         I                           bi, bj, myTime, myIter, myThid )
569            ELSEIF ( implicitViscosity ) THEN
570    #else /* INCLUDE_IMPLVERTADV_CODE */
571            IF     ( implicitViscosity ) THEN
572    #endif /* INCLUDE_IMPLVERTADV_CODE */
573  #ifdef    ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
574            idkey = iikey + 3  CADJ STORE KappaRU(:,:,:) = comlev1_bibj , key=idynkey, byte=isbyte
575  CADJ STORE gUNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE gU(:,:,:,bi,bj) = comlev1_bibj , key=idynkey, byte=isbyte
576  #endif    /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
577            CALL IMPLDIFF(            CALL IMPLDIFF(
578       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
579       I         deltaTmom, KappaRU,recip_HFacW,       I         -1, KappaRU,recip_HFacW,
580       U         gUNm1,       U         gU,
581       I         myThid )       I         myThid )
582  #ifdef    ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
583            idkey = iikey + 4  CADJ STORE KappaRV(:,:,:) = comlev1_bibj , key=idynkey, byte=isbyte
584  CADJ STORE gVNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE gV(:,:,:,bi,bj) = comlev1_bibj , key=idynkey, byte=isbyte
585  #endif    /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
586            CALL IMPLDIFF(            CALL IMPLDIFF(
587       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
588       I         deltaTmom, KappaRV,recip_HFacS,       I         -2, KappaRV,recip_HFacS,
589       U         gVNm1,       U         gV,
590       I         myThid )       I         myThid )
591            ENDIF
592    
593  #ifdef   ALLOW_OBCS  #ifdef   ALLOW_OBCS
594  C--      Apply open boundary conditions  C--      Apply open boundary conditions
595           IF (useOBCS) THEN          IF ( useOBCS .AND.(implicitViscosity.OR.momImplVertAdv) ) THEN
596             DO K=1,Nr             DO K=1,Nr
597               CALL OBCS_APPLY_UV( bi, bj, k, gUnm1, gVnm1, myThid )               CALL OBCS_APPLY_UV( bi, bj, k, gU, gV, myThid )
598             ENDDO             ENDDO
599           END IF          ENDIF
600  #endif   /* ALLOW_OBCS */  #endif   /* ALLOW_OBCS */
601    
602  #ifdef    INCLUDE_CD_CODE  #ifdef    ALLOW_CD_CODE
603            IF (implicitViscosity.AND.useCDscheme) THEN
604  #ifdef    ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
605            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  
606  #endif    /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
607            CALL IMPLDIFF(            CALL IMPLDIFF(
608       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
609       I         deltaTmom, KappaRU,recip_HFacW,       I         0, KappaRU,recip_HFacW,
610       U         vVelD,       U         vVelD,
611       I         myThid )       I         myThid )
612  #ifdef    ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
613            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  
614  #endif    /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
615            CALL IMPLDIFF(            CALL IMPLDIFF(
616       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
617       I         deltaTmom, KappaRV,recip_HFacS,       I         0, KappaRV,recip_HFacS,
618       U         uVelD,       U         uVelD,
619       I         myThid )       I         myThid )
 #endif    /* INCLUDE_CD_CODE */  
 C--     End If implicitViscosity.AND.momStepping  
620          ENDIF          ENDIF
621    #endif    /* ALLOW_CD_CODE */
622  Cjmc : add for phiHyd output <- but not working if multi tile per CPU  C--     End implicit Vertical advection & viscosity
623  c       IF ( DIFFERENT_MULTIPLE(dumpFreq,myTime+deltaTClock,myTime)  
624  c    &  .AND. buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
625  c         WRITE(suff,'(I10.10)') myIter+1  
626  c         CALL WRITE_FLD_XYZ_RL('PH.',suff,phiHyd,myIter+1,myThid)  #ifdef ALLOW_NONHYDROSTATIC
627  c       ENDIF  C--   Step forward W field in N-H algorithm
628  Cjmc(end)          IF ( nonHydrostatic ) THEN
629    #ifdef ALLOW_DEBUG
630  #ifdef ALLOW_TIMEAVE           IF ( debugLevel .GE. debLevB )
631          IF (taveFreq.GT.0.) THEN       &     CALL DEBUG_CALL('CALC_GW', myThid )
632            CALL TIMEAVE_CUMUL_1T(phiHydtave, phiHyd, Nr,  #endif
633       I                              deltaTclock, bi, bj, myThid)           CALL TIMER_START('CALC_GW          [DYNAMICS]',myThid)
634            IF (ivdc_kappa.NE.0.) THEN           CALL CALC_GW(
635              CALL TIMEAVE_CUMULATE(ConvectCountTave, ConvectCount, Nr,       I                 bi,bj, KappaRU, KappaRV,
636       I                              deltaTclock, bi, bj, myThid)       I                 myTime, myIter, myThid )
           ENDIF  
637          ENDIF          ENDIF
638  #endif /* ALLOW_TIMEAVE */          IF ( nonHydrostatic.OR.implicitIntGravWave )
639         &   CALL TIMESTEP_WVEL( bi,bj, myTime, myIter, myThid )
640            IF ( nonHydrostatic )
641         &   CALL TIMER_STOP ('CALC_GW          [DYNAMICS]',myThid)
642    #endif
643    
644    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
645    
646    C-    end of bi,bj loops
647         ENDDO         ENDDO
648        ENDDO        ENDDO
649    
650  #ifndef EXCLUDE_DEBUGMODE  #ifdef ALLOW_OBCS
651        If (debugMode) THEN        IF (useOBCS) THEN
652           CALL OBCS_PRESCRIBE_EXCHANGES(myThid)
653          ENDIF
654    #endif
655    
656    Cml(
657    C     In order to compare the variance of phiHydLow of a p/z-coordinate
658    C     run with etaH of a z/p-coordinate run the drift of phiHydLow
659    C     has to be removed by something like the following subroutine:
660    C      CALL REMOVE_MEAN_RL( 1, phiHydLow, maskInC, maskInC, rA, drF,
661    C     &                     'phiHydLow', myTime, myThid )
662    Cml)
663    
664    #ifdef ALLOW_DIAGNOSTICS
665          IF ( useDiagnostics ) THEN
666    
667           CALL DIAGNOSTICS_FILL(totPhihyd,'PHIHYD  ',0,Nr,0,1,1,myThid)
668           CALL DIAGNOSTICS_FILL(phiHydLow,'PHIBOT  ',0, 1,0,1,1,myThid)
669    
670           tmpFac = 1. _d 0
671           CALL DIAGNOSTICS_SCALE_FILL(totPhihyd,tmpFac,2,
672         &                                 'PHIHYDSQ',0,Nr,0,1,1,myThid)
673    
674           CALL DIAGNOSTICS_SCALE_FILL(phiHydLow,tmpFac,2,
675         &                                 'PHIBOTSQ',0, 1,0,1,1,myThid)
676    
677          ENDIF
678    #endif /* ALLOW_DIAGNOSTICS */
679    
680    #ifdef ALLOW_DEBUG
681          If ( debugLevel .GE. debLevB ) THEN
682         CALL DEBUG_STATS_RL(1,EtaN,'EtaN (DYNAMICS)',myThid)         CALL DEBUG_STATS_RL(1,EtaN,'EtaN (DYNAMICS)',myThid)
683           CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (DYNAMICS)',myThid)
684         CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (DYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (DYNAMICS)',myThid)
685         CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (DYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (DYNAMICS)',myThid)
686         CALL DEBUG_STATS_RL(Nr,theta,'Theta (DYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,theta,'Theta (DYNAMICS)',myThid)
687         CALL DEBUG_STATS_RL(Nr,salt,'Salt (DYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,salt,'Salt (DYNAMICS)',myThid)
688         CALL DEBUG_STATS_RL(Nr,Gu,'Gu (DYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,gU,'Gu (DYNAMICS)',myThid)
689         CALL DEBUG_STATS_RL(Nr,Gv,'Gv (DYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,gV,'Gv (DYNAMICS)',myThid)
690         CALL DEBUG_STATS_RL(Nr,Gt,'Gt (DYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,gT,'Gt (DYNAMICS)',myThid)
691         CALL DEBUG_STATS_RL(Nr,Gs,'Gs (DYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,gS,'Gs (DYNAMICS)',myThid)
692         CALL DEBUG_STATS_RL(Nr,GuNm1,'GuNm1 (DYNAMICS)',myThid)  #ifndef ALLOW_ADAMSBASHFORTH_3
693         CALL DEBUG_STATS_RL(Nr,GvNm1,'GvNm1 (DYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,guNm1,'GuNm1 (DYNAMICS)',myThid)
694         CALL DEBUG_STATS_RL(Nr,GtNm1,'GtNm1 (DYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,gvNm1,'GvNm1 (DYNAMICS)',myThid)
695         CALL DEBUG_STATS_RL(Nr,GsNm1,'GsNm1 (DYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,gtNm1,'GtNm1 (DYNAMICS)',myThid)
696           CALL DEBUG_STATS_RL(Nr,gsNm1,'GsNm1 (DYNAMICS)',myThid)
697    #endif
698          ENDIF
699    #endif
700    
701    #ifdef DYNAMICS_GUGV_EXCH_CHECK
702    C- jmc: For safety checking only: This Exchange here should not change
703    C       the solution. If solution changes, it means something is wrong,
704    C       but it does not mean that it is less wrong with this exchange.
705          IF ( debugLevel .GT. debLevB ) THEN
706           CALL EXCH_UV_XYZ_RL(gU,gV,.TRUE.,myThid)
707        ENDIF        ENDIF
708  #endif  #endif
709    
710    #ifdef ALLOW_DEBUG
711          IF ( debugLevel .GE. debLevB )
712         &   CALL DEBUG_LEAVE( 'DYNAMICS', myThid )
713    #endif
714    
715        RETURN        RETURN
716        END        END
Legend:
Removed from v.1.70  
changed lines
  Added in v.1.146
  ViewVC Help
Powered by ViewVC 1.1.22