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

Legend:
Removed from v.1.69  
changed lines
  Added in v.1.161

  ViewVC Help
Powered by ViewVC 1.1.22