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