/[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.67 by heimbach, Mon May 14 21:46:17 2001 UTC revision 1.168 by jmc, Fri Jan 3 16:19:04 2014 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_MOM_COMMON
7    # include "MOM_COMMON_OPTIONS.h"
8    #endif
9    #ifdef ALLOW_OBCS
10    # include "OBCS_OPTIONS.h"
11    #endif
12    
13    #undef DYNAMICS_GUGV_EXCH_CHECK
14    
15    CBOP
16    C     !ROUTINE: DYNAMICS
17    C     !INTERFACE:
18        SUBROUTINE DYNAMICS(myTime, myIter, myThid)        SUBROUTINE DYNAMICS(myTime, myIter, myThid)
19  C     /==========================================================\  C     !DESCRIPTION: \bv
20  C     | SUBROUTINE DYNAMICS                                      |  C     *==========================================================*
21  C     | o Controlling routine for the explicit part of the model |  C     | SUBROUTINE DYNAMICS
22  C     |   dynamics.                                              |  C     | o Controlling routine for the explicit part of the model
23  C     |==========================================================|  C     |   dynamics.
24  C     | This routine evaluates the "dynamics" terms for each     |  C     *==========================================================*
25  C     | block of ocean in turn. Because the blocks of ocean have |  C     | This routine evaluates the "dynamics" terms for each
26  C     | overlap regions they are independent of one another.     |  C     | block of ocean in turn. Because the blocks of ocean have
27  C     | If terms involving lateral integrals are needed in this  |  C     | overlap regions they are independent of one another.
28  C     | routine care will be needed. Similarly finite-difference |  C     | If terms involving lateral integrals are needed in this
29  C     | operations with stencils wider than the overlap region   |  C     | routine care will be needed. Similarly finite-difference
30  C     | require special consideration.                           |  C     | operations with stencils wider than the overlap region
31  C     | Notes                                                    |  C     | require special consideration.
32  C     | =====                                                    |  C     | The algorithm...
33  C     | C*P* comments indicating place holders for which code is |  C     |
34  C     |      presently being developed.                          |  C     | "Correction Step"
35  C     \==========================================================/  C     | =================
36    C     | Here we update the horizontal velocities with the surface
37    C     | pressure such that the resulting flow is either consistent
38    C     | with the free-surface evolution or the rigid-lid:
39    C     |   U[n] = U* + dt x d/dx P
40    C     |   V[n] = V* + dt x d/dy P
41    C     |   W[n] = W* + dt x d/dz P  (NH mode)
42    C     |
43    C     | "Calculation of Gs"
44    C     | ===================
45    C     | This is where all the accelerations and tendencies (ie.
46    C     | physics, parameterizations etc...) are calculated
47    C     |   rho = rho ( theta[n], salt[n] )
48    C     |   b   = b(rho, theta)
49    C     |   K31 = K31 ( rho )
50    C     |   Gu[n] = Gu( u[n], v[n], wVel, b, ... )
51    C     |   Gv[n] = Gv( u[n], v[n], wVel, b, ... )
52    C     |   Gt[n] = Gt( theta[n], u[n], v[n], wVel, K31, ... )
53    C     |   Gs[n] = Gs( salt[n], u[n], v[n], wVel, K31, ... )
54    C     |
55    C     | "Time-stepping" or "Prediction"
56    C     | ================================
57    C     | The models variables are stepped forward with the appropriate
58    C     | time-stepping scheme (currently we use Adams-Bashforth II)
59    C     | - For momentum, the result is always *only* a "prediction"
60    C     | in that the flow may be divergent and will be "corrected"
61    C     | later with a surface pressure gradient.
62    C     | - Normally for tracers the result is the new field at time
63    C     | level [n+1} *BUT* in the case of implicit diffusion the result
64    C     | is also *only* a prediction.
65    C     | - We denote "predictors" with an asterisk (*).
66    C     |   U* = U[n] + dt x ( 3/2 Gu[n] - 1/2 Gu[n-1] )
67    C     |   V* = V[n] + dt x ( 3/2 Gv[n] - 1/2 Gv[n-1] )
68    C     |   theta[n+1] = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
69    C     |   salt[n+1] = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
70    C     | With implicit diffusion:
71    C     |   theta* = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
72    C     |   salt* = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
73    C     |   (1 + dt * K * d_zz) theta[n] = theta*
74    C     |   (1 + dt * K * d_zz) salt[n] = salt*
75    C     |
76    C     *==========================================================*
77    C     \ev
78    C     !USES:
79        IMPLICIT NONE        IMPLICIT NONE
   
80  C     == Global variables ===  C     == Global variables ===
81  #include "SIZE.h"  #include "SIZE.h"
82  #include "EEPARAMS.h"  #include "EEPARAMS.h"
83  #include "PARAMS.h"  #include "PARAMS.h"
 #include "DYNVARS.h"  
84  #include "GRID.h"  #include "GRID.h"
85    #include "DYNVARS.h"
86    #ifdef ALLOW_MOM_COMMON
87    # include "MOM_VISC.h"
88    #endif
89    #ifdef ALLOW_CD_CODE
90    # include "CD_CODE_VARS.h"
91    #endif
92  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
93  # include "tamc.h"  # include "tamc.h"
94  # include "tamc_keys.h"  # include "tamc_keys.h"
95  # include "FFIELDS.h"  # include "FFIELDS.h"
96    # include "EOS.h"
97  # ifdef ALLOW_KPP  # ifdef ALLOW_KPP
98  #  include "KPP.h"  #  include "KPP.h"
99  # endif  # endif
100  # ifdef ALLOW_GMREDI  # ifdef ALLOW_PTRACERS
101  #  include "GMREDI.h"  #  include "PTRACERS_SIZE.h"
102    #  include "PTRACERS_FIELDS.h"
103    # endif
104    # ifdef ALLOW_OBCS
105    #  include "OBCS_PARAMS.h"
106    #  include "OBCS_FIELDS.h"
107    #  ifdef ALLOW_PTRACERS
108    #   include "OBCS_PTRACERS.h"
109    #  endif
110    # endif
111    # ifdef ALLOW_MOM_FLUXFORM
112    #  include "MOM_FLUXFORM.h"
113  # endif  # endif
114  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
115    
116  #ifdef ALLOW_TIMEAVE  C     !CALLING SEQUENCE:
117  #include "TIMEAVE_STATV.h"  C     DYNAMICS()
118  #endif  C      |
119    C      |-- CALC_EP_FORCING
120    C      |
121    C      |-- CALC_GRAD_PHI_SURF
122    C      |
123    C      |-- CALC_VISCOSITY
124    C      |
125    C      |-- MOM_CALC_3D_STRAIN
126    C      |
127    C      |-- CALC_EDDY_STRESS
128    C      |
129    C      |-- CALC_PHI_HYD
130    C      |
131    C      |-- MOM_FLUXFORM
132    C      |
133    C      |-- MOM_VECINV
134    C      |
135    C      |-- MOM_CALC_SMAG_3D
136    C      |-- MOM_UV_SMAG_3D
137    C      |
138    C      |-- TIMESTEP
139    C      |
140    C      |-- MOM_U_IMPLICIT_R
141    C      |-- MOM_V_IMPLICIT_R
142    C      |
143    C      |-- IMPLDIFF
144    C      |
145    C      |-- OBCS_APPLY_UV
146    C      |
147    C      |-- CALC_GW
148    C      |
149    C      |-- DIAGNOSTICS_FILL
150    C      |-- DEBUG_STATS_RL
151    
152    C     !INPUT/OUTPUT PARAMETERS:
153  C     == Routine arguments ==  C     == Routine arguments ==
154  C     myTime - Current time in simulation  C     myTime :: Current time in simulation
155  C     myIter - Current iteration number in simulation  C     myIter :: Current iteration number in simulation
156  C     myThid - Thread number for this instance of the routine.  C     myThid :: Thread number for this instance of the routine.
157        _RL myTime        _RL myTime
158        INTEGER myIter        INTEGER myIter
159        INTEGER myThid        INTEGER myThid
160    
161    C     !FUNCTIONS:
162    #ifdef ALLOW_DIAGNOSTICS
163          LOGICAL  DIAGNOSTICS_IS_ON
164          EXTERNAL DIAGNOSTICS_IS_ON
165    #endif
166    
167    C     !LOCAL VARIABLES:
168  C     == Local variables  C     == Local variables
169  C     xA, yA                 - Per block temporaries holding face areas  C     fVer[UV]               o fVer: Vertical flux term - note fVer
170  C     uTrans, vTrans, rTrans - Per block temporaries holding flow  C                                    is "pipelined" in the vertical
171  C                              transport  C                                    so we need an fVer for each
172  C                              o uTrans: Zonal transport  C                                    variable.
173  C                              o vTrans: Meridional transport  C     phiHydC    :: hydrostatic potential anomaly at cell center
174  C                              o rTrans: Vertical transport  C                   In z coords phiHyd is the hydrostatic potential
175  C     maskC,maskUp             o maskC: land/water mask for tracer cells  C                      (=pressure/rho0) anomaly
176  C                              o maskUp: land/water mask for W points  C                   In p coords phiHyd is the geopotential height anomaly.
177  C     fVer[STUV]               o fVer: Vertical flux term - note fVer  C     phiHydF    :: hydrostatic potential anomaly at middle between 2 centers
178  C                                      is "pipelined" in the vertical  C     dPhiHydX,Y :: Gradient (X & Y directions) of hydrostatic potential anom.
179  C                                      so we need an fVer for each  C     phiSurfX,  ::  gradient of Surface potential (Pressure/rho, ocean)
180  C                                      variable.  C     phiSurfY             or geopotential (atmos) in X and Y direction
181  C     rhoK, rhoKM1   - Density at current level, and level above  C     guDissip   :: dissipation tendency (all explicit terms), u component
182  C     phiHyd         - Hydrostatic part of the potential phiHydi.  C     gvDissip   :: dissipation tendency (all explicit terms), v component
183  C                      In z coords phiHydiHyd is the hydrostatic  C     KappaRU    :: vertical viscosity for velocity U-component
184  C                      Potential (=pressure/rho0) anomaly  C     KappaRV    :: vertical viscosity for velocity V-component
185  C                      In p coords phiHydiHyd is the geopotential  C     iMin, iMax :: Ranges and sub-block indices on which calculations
186  C                      surface height anomaly.  C     jMin, jMax    are applied.
187  C     phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean)  C     bi, bj     :: tile indices
188  C     phiSurfY             or geopotentiel (atmos) in X and Y direction  C     k          :: current level index
189  C     KappaRT,       - Total diffusion in vertical for T and S.  C     km1, kp1   :: index of level above (k-1) and below (k+1)
190  C     KappaRS          (background + spatially varying, isopycnal term).  C     kUp, kDown :: Index for interface above and below. kUp and kDown are
191  C     iMin, iMax     - Ranges and sub-block indices on which calculations  C                   are switched with k to be the appropriate index into fVerU,V
 C     jMin, jMax       are applied.  
 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.  
       _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 maskC   (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)  
192        _RL fVerU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
193        _RL fVerV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
194        _RL phiHyd  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL phiHydF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
195        _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiHydC (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
196        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL dPhiHydX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
197          _RL dPhiHydY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
198        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
199        _RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
200        _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL guDissip(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
201        _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL gvDissip(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
202        _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
203        _RL KappaRV (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
204        _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  #ifdef ALLOW_SMAG_3D
205        _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  C     str11       :: strain component Vxx @ grid-cell center
206        _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  C     str22       :: strain component Vyy @ grid-cell center
207    C     str33       :: strain component Vzz @ grid-cell center
208  C This is currently used by IVDC and Diagnostics  C     str12       :: strain component Vxy @ grid-cell corner
209        _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  C     str13       :: strain component Vxz @ above uVel
210    C     str23       :: strain component Vyz @ above vVel
211    C     viscAh3d_00 :: Smagorinsky viscosity @ grid-cell center
212    C     viscAh3d_12 :: Smagorinsky viscosity @ grid-cell corner
213    C     viscAh3d_13 :: Smagorinsky viscosity @ above uVel
214    C     viscAh3d_23 :: Smagorinsky viscosity @ above vVel
215    C     addDissU    :: zonal momentum tendency from 3-D Smag. viscosity
216    C     addDissV    :: merid momentum tendency from 3-D Smag. viscosity
217          _RL str11(1-OLx:sNx+OLx,1-OLy:sNy+OLy, Nr )
218          _RL str22(1-OLx:sNx+OLx,1-OLy:sNy+OLy, Nr )
219          _RL str33(1-OLx:sNx+OLx,1-OLy:sNy+OLy, Nr )
220          _RL str12(1-OLx:sNx+OLx,1-OLy:sNy+OLy, Nr )
221          _RL str13(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)
222          _RL str23(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)
223          _RL viscAh3d_00(1-OLx:sNx+OLx,1-OLy:sNy+OLy, Nr )
224          _RL viscAh3d_12(1-OLx:sNx+OLx,1-OLy:sNy+OLy, Nr )
225          _RL viscAh3d_13(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)
226          _RL viscAh3d_23(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)
227          _RL addDissU(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
228          _RL addDissV(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
229    #elif ( defined ALLOW_NONHYDROSTATIC )
230          _RL str13(1), str23(1), str33(1)
231          _RL viscAh3d_00(1), viscAh3d_13(1), viscAh3d_23(1)
232    #endif
233    
       INTEGER iMin, iMax  
       INTEGER jMin, jMax  
234        INTEGER bi, bj        INTEGER bi, bj
235        INTEGER i, j        INTEGER i, j
236        INTEGER k, km1, kup, kDown        INTEGER k, km1, kp1, kUp, kDown
237          INTEGER iMin, iMax
238          INTEGER jMin, jMax
239          PARAMETER( iMin = 0 , iMax = sNx+1 )
240          PARAMETER( jMin = 0 , jMax = sNy+1 )
241    
242    #ifdef ALLOW_DIAGNOSTICS
243          LOGICAL dPhiHydDiagIsOn
244          _RL tmpFac
245    #endif /* ALLOW_DIAGNOSTICS */
246    
 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)  
   
247  C---    The algorithm...  C---    The algorithm...
248  C  C
249  C       "Correction Step"  C       "Correction Step"
# Line 165  C         salt* = salt[n] + dt x ( 3/2 G Line 287  C         salt* = salt[n] + dt x ( 3/2 G
287  C         (1 + dt * K * d_zz) theta[n] = theta*  C         (1 + dt * K * d_zz) theta[n] = theta*
288  C         (1 + dt * K * d_zz) salt[n] = salt*  C         (1 + dt * K * d_zz) salt[n] = salt*
289  C---  C---
290    CEOP
291    
292  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_DEBUG
293  C--   dummy statement to end declaration part        IF (debugMode) CALL DEBUG_ENTER( 'DYNAMICS', myThid )
294        ikey = 1  #endif
 #endif /* ALLOW_AUTODIFF_TAMC */  
295    
296  C--   Set up work arrays with valid (i.e. not NaN) values  #ifdef ALLOW_DIAGNOSTICS
297  C     These inital values do not alter the numerical results. They        dPhiHydDiagIsOn = .FALSE.
298  C     just ensure that all memory references are to valid floating        IF ( useDiagnostics )
299  C     point numbers. This prevents spurious hardware signals due to       &  dPhiHydDiagIsOn = DIAGNOSTICS_IS_ON( 'Um_dPHdx', myThid )
300  C     uninitialised but inert locations.       &               .OR. DIAGNOSTICS_IS_ON( 'Vm_dPHdy', myThid )
301        DO j=1-OLy,sNy+OLy  #endif
302         DO i=1-OLx,sNx+OLx  
303          xA(i,j)      = 0. _d 0  C-- Call to routine for calculation of Eliassen-Palm-flux-forced
304          yA(i,j)      = 0. _d 0  C    U-tendency, if desired:
305          uTrans(i,j)  = 0. _d 0  #ifdef INCLUDE_EP_FORCING_CODE
306          vTrans(i,j)  = 0. _d 0        CALL CALC_EP_FORCING(myThid)
307          DO k=1,Nr  #endif
          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  
         maskC  (i,j) = 0. _d 0  
         phiSurfX(i,j) = 0. _d 0  
         phiSurfY(i,j) = 0. _d 0  
        ENDDO  
       ENDDO  
308    
309    #ifdef ALLOW_AUTODIFF_MONITOR_DIAG
310          CALL DUMMY_IN_DYNAMICS( myTime, myIter, myThid )
311    #endif
312    
313  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
314  C--   HPF directive to help TAMC  C--   HPF directive to help TAMC
# Line 208  CHPF$ INDEPENDENT Line 319  CHPF$ INDEPENDENT
319    
320  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
321  C--    HPF directive to help TAMC  C--    HPF directive to help TAMC
322  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS,fVerU,fVerV  CHPF$  INDEPENDENT, NEW (fVerU,fVerV
323  CHPF$&                  ,phiHyd,utrans,vtrans,maskc,xA,yA  CHPF$&                  ,phiHydF
324  CHPF$&                  ,KappaRT,KappaRS,KappaRU,KappaRV  CHPF$&                  ,KappaRU,KappaRV
325  CHPF$&                  )  CHPF$&                  )
326  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
327    
# Line 219  CHPF$&                  ) Line 330  CHPF$&                  )
330  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
331            act1 = bi - myBxLo(myThid)            act1 = bi - myBxLo(myThid)
332            max1 = myBxHi(myThid) - myBxLo(myThid) + 1            max1 = myBxHi(myThid) - myBxLo(myThid) + 1
   
333            act2 = bj - myByLo(myThid)            act2 = bj - myByLo(myThid)
334            max2 = myByHi(myThid) - myByLo(myThid) + 1            max2 = myByHi(myThid) - myByLo(myThid) + 1
   
335            act3 = myThid - 1            act3 = myThid - 1
336            max3 = nTx*nTy            max3 = nTx*nTy
   
337            act4 = ikey_dynamics - 1            act4 = ikey_dynamics - 1
338              idynkey = (act1 + 1) + act2*max1
           ikey = (act1 + 1) + act2*max1  
339       &                      + act3*max1*max2       &                      + act3*max1*max2
340       &                      + act4*max1*max2*max3       &                      + act4*max1*max2*max3
341  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
342    
343  C--     Set up work arrays that need valid initial values  C--   Set up work arrays with valid (i.e. not NaN) values
344          DO j=1-OLy,sNy+OLy  C     These initial values do not alter the numerical results. They
345           DO i=1-OLx,sNx+OLx  C     just ensure that all memory references are to valid floating
346            rTrans(i,j)   = 0. _d 0  C     point numbers. This prevents spurious hardware signals due to
347            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  
348    
349    #ifdef ALLOW_AUTODIFF
350          DO k=1,Nr          DO k=1,Nr
351           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
352            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
353  C This is currently also used by IVDC and Diagnostics  c--   need some re-initialisation here to break dependencies
354             ConvectCount(i,j,k) = 0.             gU(i,j,k,bi,bj) = 0. _d 0
355             KappaRT(i,j,k) = 0. _d 0             gV(i,j,k,bi,bj) = 0. _d 0
            KappaRS(i,j,k) = 0. _d 0  
356            ENDDO            ENDDO
357           ENDDO           ENDDO
358          ENDDO          ENDDO
359    #endif /* ALLOW_AUTODIFF */
360          iMin = 1-OLx+1          DO j=1-OLy,sNy+OLy
361          iMax = sNx+OLx           DO i=1-OLx,sNx+OLx
362          jMin = 1-OLy+1            fVerU  (i,j,1) = 0. _d 0
363          jMax = sNy+OLy            fVerU  (i,j,2) = 0. _d 0
364              fVerV  (i,j,1) = 0. _d 0
365              fVerV  (i,j,2) = 0. _d 0
366  #ifdef ALLOW_AUTODIFF_TAMC            phiHydF (i,j)  = 0. _d 0
367  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte            phiHydC (i,j)  = 0. _d 0
368  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  #ifndef INCLUDE_PHIHYD_CALCULATION_CODE
369  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte            dPhiHydX(i,j)  = 0. _d 0
370  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte            dPhiHydY(i,j)  = 0. _d 0
371  #endif /* ALLOW_AUTODIFF_TAMC */  #endif
372              phiSurfX(i,j)  = 0. _d 0
373  C--     Start of diagnostic loop            phiSurfY(i,j)  = 0. _d 0
374          DO k=Nr,1,-1            guDissip(i,j)  = 0. _d 0
375              gvDissip(i,j)  = 0. _d 0
376  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF
377  C? Patrick, is this formula correct now that we change the loop range?            phiHydLow(i,j,bi,bj) = 0. _d 0
378  C? Do we still need this?  # if (defined NONLIN_FRSURF) && (defined ALLOW_MOM_FLUXFORM)
379  cph kkey formula corrected.  #  ifndef DISABLE_RSTAR_CODE
380  cph Needed for rhok, rhokm1, in the case useGMREDI.            dWtransC(i,j,bi,bj) = 0. _d 0
381           kkey = (ikey-1)*Nr + k            dWtransU(i,j,bi,bj) = 0. _d 0
382  CADJ STORE rhokm1(:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte            dWtransV(i,j,bi,bj) = 0. _d 0
383  CADJ STORE rhok  (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte  #  endif
384  #endif /* ALLOW_AUTODIFF_TAMC */  # endif
385    #endif /* ALLOW_AUTODIFF */
386  C--       Integrate continuity vertically for vertical velocity           ENDDO
           CALL INTEGRATE_FOR_W(  
      I                         bi, bj, k, uVel, vVel,  
      O                         wVel,  
      I                         myThid )  
   
 #ifdef    ALLOW_OBCS  
 #ifdef    ALLOW_NONHYDROSTATIC  
 C--       Apply OBC to W if in N-H mode  
           IF (useOBCS.AND.nonHydrostatic) THEN  
             CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )  
           ENDIF  
 #endif    /* ALLOW_NONHYDROSTATIC */  
 #endif    /* ALLOW_OBCS */  
   
 C--       Calculate gradients of potential density for isoneutral  
 C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)  
 c         IF ( k.GT.1 .AND. (useGMRedi.OR.ivdc_kappa.NE.0.) ) THEN  
           IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.) ) THEN  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  
 CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
             CALL FIND_RHO(  
      I        bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,  
      I        theta, salt,  
      O        rhoK,  
      I        myThid )  
             IF (k.GT.1) THEN  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  
 CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
              CALL FIND_RHO(  
      I        bi, bj, iMin, iMax, jMin, jMax, k-1, k, eosType,  
      I        theta, salt,  
      O        rhoKm1,  
      I        myThid )  
             ENDIF  
             CALL GRAD_SIGMA(  
      I             bi, bj, iMin, iMax, jMin, jMax, k,  
      I             rhoK, rhoKm1, rhoK,  
      O             sigmaX, sigmaY, sigmaR,  
      I             myThid )  
           ENDIF  
   
 C--       Implicit Vertical Diffusion for Convection  
 c ==> should use sigmaR !!!  
           IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN  
             CALL CALC_IVDC(  
      I        bi, bj, iMin, iMax, jMin, jMax, k,  
      I        rhoKm1, rhoK,  
      U        ConvectCount, KappaRT, KappaRS,  
      I        myTime, myIter, myThid)  
           ENDIF  
   
 C--     end of diagnostic k loop (Nr:1)  
387          ENDDO          ENDDO
388    
389    C--     Start computation of dynamics
390    
391  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
392  cph avoids recomputation of integrate_for_w  CADJ STORE wVel (:,:,:,bi,bj) =
393  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ &     comlev1_bibj, key=idynkey, byte=isbyte
394  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
395    
396  #ifdef  ALLOW_OBCS  C--     Explicit part of the Surface Potential Gradient (add in TIMESTEP)
397  C--     Calculate future values on open boundaries  C       (note: this loop will be replaced by CALL CALC_GRAD_ETA)
398          IF (useOBCS) THEN          IF (implicSurfPress.NE.1.) THEN
399            CALL OBCS_CALC( bi, bj, myTime+deltaT,            CALL CALC_GRAD_PHI_SURF(
400       I            uVel, vVel, wVel, theta, salt,       I         bi,bj,iMin,iMax,jMin,jMax,
401       I            myThid )       I         etaN,
402         O         phiSurfX,phiSurfY,
403         I         myThid )
404          ENDIF          ENDIF
 #endif  /* ALLOW_OBCS */  
405    
 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 )  
406  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
407  cph needed for KPP  CADJ STORE uVel (:,:,:,bi,bj) = comlev1_bibj, key=idynkey, byte=isbyte
408  CADJ STORE surfacetendencyU(:,:,bi,bj)  CADJ STORE vVel (:,:,:,bi,bj) = comlev1_bibj, key=idynkey, byte=isbyte
409  CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  #ifdef ALLOW_KPP
410  CADJ STORE surfacetendencyV(:,:,bi,bj)  CADJ STORE KPPviscAz (:,:,:,bi,bj)
411  CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  CADJ &                 = comlev1_bibj, key=idynkey, byte=isbyte
412  CADJ STORE surfacetendencyS(:,:,bi,bj)  #endif /* ALLOW_KPP */
 CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE surfacetendencyT(:,:,bi,bj)  
 CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  
413  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
414    
415  #ifdef  ALLOW_GMREDI  #ifndef ALLOW_AUTODIFF
416            IF ( .NOT.momViscosity ) THEN
417  #ifdef ALLOW_AUTODIFF_TAMC  #endif
 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  
418            DO k=1,Nr            DO k=1,Nr
419              CALL GMREDI_CALC_TENSOR(             DO j=1-OLy,sNy+OLy
420       I             bi, bj, iMin, iMax, jMin, jMax, k,              DO i=1-OLx,sNx+OLx
421       I             sigmaX, sigmaY, sigmaR,               KappaRU(i,j,k) = 0. _d 0
422       I             myThid )               KappaRV(i,j,k) = 0. _d 0
423            ENDDO              ENDDO
424  #ifdef ALLOW_AUTODIFF_TAMC             ENDDO
         ELSE  
           DO k=1, Nr  
             CALL GMREDI_CALC_TENSOR_DUMMY(  
      I             bi, bj, iMin, iMax, jMin, jMax, k,  
      I             sigmaX, sigmaY, sigmaR,  
      I             myThid )  
425            ENDDO            ENDDO
426  #endif /* ALLOW_AUTODIFF_TAMC */  #ifndef ALLOW_AUTODIFF
427          ENDIF          ENDIF
428    #endif
429  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef INCLUDE_CALC_DIFFUSIVITY_CALL
430  CADJ STORE Kwx(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte  C--     Calculate the total vertical viscosity
431  CADJ STORE Kwy(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte          IF ( momViscosity ) THEN
432  CADJ STORE Kwz(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte            CALL CALC_VISCOSITY(
433  #endif /* ALLOW_AUTODIFF_TAMC */       I            bi,bj, iMin,iMax,jMin,jMax,
434         O            KappaRU, KappaRV,
435  #endif  /* ALLOW_GMREDI */       I            myThid )
   
 #ifdef  ALLOW_KPP  
 C--     Compute KPP mixing coefficients  
         IF (useKPP) THEN  
           CALL KPP_CALC(  
      I                  bi, bj, myTime, myThid )  
 #ifdef ALLOW_AUTODIFF_TAMC  
         ELSE  
           CALL KPP_CALC_DUMMY(  
      I                  bi, bj, myTime, myThid )  
 #endif /* ALLOW_AUTODIFF_TAMC */  
436          ENDIF          ENDIF
437    #endif /* INCLUDE_CALC_DIFFUSIVITY_CALL */
438    
439  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_SMAG_3D
440  CADJ STORE KPPghat   (:,:,:,bi,bj)          IF ( useSmag3D ) THEN
441  CADJ &   , KPPviscAz (:,:,:,bi,bj)            CALL MOM_CALC_3D_STRAIN(
442  CADJ &   , KPPdiffKzT(:,:,:,bi,bj)       O         str11, str22, str33, str12, str13, str23,
443  CADJ &   , KPPdiffKzS(:,:,:,bi,bj)       I         bi, bj, myThid )
 CADJ &   , KPPfrac   (:,:  ,bi,bj)  
 CADJ &                 = comlev1_bibj, key=ikey, byte=isbyte  
 #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)  
444          ENDIF          ENDIF
445  #endif /* ALLOW_AIM */  #endif /* ALLOW_SMAG_3D */
   
   
 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,km1,kup,kDown,  
      O        xA,yA,uTrans,vTrans,rTrans,maskC,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  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL  
 C--      Calculate the total vertical diffusivity  
          CALL CALC_DIFFUSIVITY(  
      I        bi,bj,iMin,iMax,jMin,jMax,k,  
      I        maskC,maskup,  
      O        KappaRT,KappaRS,KappaRU,KappaRV,  
      I        myThid)  
 #endif  
   
 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,maskC,  
      I         KappaRT,  
      U         fVerT,  
      I         myTime, myThid)  
            CALL TIMESTEP_TRACER(  
      I         bi,bj,iMin,iMax,jMin,jMax,k,  
      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,maskC,  
      I         KappaRS,  
      U         fVerS,  
      I         myTime, myThid)  
            CALL TIMESTEP_TRACER(  
      I         bi,bj,iMin,iMax,jMin,jMax,k,  
      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)  
         ENDDO  
   
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 C? Patrick? What about this one?  
 cph Keys iikey and idkey don't seem to be needed  
 cph since storing occurs on different tape for each  
 cph impldiff call anyways.  
 cph Thus, common block comlev1_impl isn't needed either.  
 cph Storing below needed in the case useGMREDI.  
         iikey = (ikey-1)*maximpl  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 C--     Implicit diffusion  
         IF (implicitDiffusion) THEN  
446    
          IF (tempStepping) THEN  
447  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
448              idkey = iikey + 1  CADJ STORE KappaRU(:,:,:)
449  CADJ STORE gTNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ &     = comlev1_bibj, key=idynkey, byte=isbyte
450  #endif /* ALLOW_AUTODIFF_TAMC */  CADJ STORE KappaRV(:,:,:)
451              CALL IMPLDIFF(  CADJ &     = comlev1_bibj, key=idynkey, byte=isbyte
452       I         bi, bj, iMin, iMax, jMin, jMax,  #endif /* ALLOW_AUTODIFF_TAMC */
453       I         deltaTtracer, KappaRT, recip_HFacC,  
454       U         gTNm1,  #ifdef ALLOW_OBCS
455       I         myThid )  C--   For Stevens boundary conditions velocities need to be extrapolated
456           ENDIF  C     (copied) to a narrow strip outside the domain
457            IF ( useOBCS ) THEN
458           IF (saltStepping) THEN            CALL OBCS_COPY_UV_N(
459  #ifdef ALLOW_AUTODIFF_TAMC       U         uVel(1-OLx,1-OLy,1,bi,bj),
460           idkey = iikey + 2       U         vVel(1-OLx,1-OLy,1,bi,bj),
461  CADJ STORE gSNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte       I         Nr, bi, bj, myThid )
 #endif /* ALLOW_AUTODIFF_TAMC */  
             CALL IMPLDIFF(  
      I         bi, bj, iMin, iMax, jMin, jMax,  
      I         deltaTtracer, KappaRS, recip_HFacC,  
      U         gSNm1,  
      I         myThid )  
          ENDIF  
   
 #ifdef   ALLOW_OBCS  
 C--      Apply open boundary conditions  
          IF (useOBCS) THEN  
            DO K=1,Nr  
              CALL OBCS_APPLY_TS( bi, bj, k, gTnm1, gSnm1, myThid )  
            ENDDO  
          END IF  
 #endif   /* ALLOW_OBCS */  
   
 C--     End If implicitDiffusion  
462          ENDIF          ENDIF
463    #endif /* ALLOW_OBCS */
464    
465  C--     Start computation of dynamics  #ifdef ALLOW_EDDYPSI
466          iMin = 1-OLx+2          CALL CALC_EDDY_STRESS(bi,bj,myThid)
467          iMax = sNx+OLx-1  #endif
         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  
468    
469  C--     Start of dynamics loop  C--     Start of dynamics loop
470          DO k=1,Nr          DO k=1,Nr
# Line 616  C--       kup    Cycles through 1,2 to p Line 474  C--       kup    Cycles through 1,2 to p
474  C--       kDown  Cycles through 2,1 to point to current layer  C--       kDown  Cycles through 2,1 to point to current layer
475    
476            km1  = MAX(1,k-1)            km1  = MAX(1,k-1)
477              kp1  = MIN(k+1,Nr)
478            kup  = 1+MOD(k+1,2)            kup  = 1+MOD(k+1,2)
479            kDown= 1+MOD(k,2)            kDown= 1+MOD(k,2)
480    
481  C--      Integrate hydrostatic balance for phiHyd with BC of  #ifdef ALLOW_AUTODIFF_TAMC
482  C        phiHyd(z=0)=0           kkey = (idynkey-1)*Nr + k
483  C        distinguishe between Stagger and Non Stagger time stepping  CADJ STORE totPhiHyd (:,:,k,bi,bj)
484           IF (staggerTimeStep) THEN  CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
485    CADJ STORE theta (:,:,k,bi,bj)
486    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
487    CADJ STORE salt  (:,:,k,bi,bj)
488    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
489    CADJ STORE gT(:,:,k,bi,bj)
490    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
491    CADJ STORE gS(:,:,k,bi,bj)
492    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
493    # ifdef NONLIN_FRSURF
494    cph-test
495    CADJ STORE  phiHydC (:,:)
496    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
497    CADJ STORE  phiHydF (:,:)
498    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
499    CADJ STORE gU(:,:,k,bi,bj)
500    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
501    CADJ STORE gV(:,:,k,bi,bj)
502    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
503    #  ifndef ALLOW_ADAMSBASHFORTH_3
504    CADJ STORE guNm1(:,:,k,bi,bj)
505    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
506    CADJ STORE gvNm1(:,:,k,bi,bj)
507    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
508    #  else
509    CADJ STORE guNm(:,:,k,bi,bj,1)
510    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
511    CADJ STORE guNm(:,:,k,bi,bj,2)
512    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
513    CADJ STORE gvNm(:,:,k,bi,bj,1)
514    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
515    CADJ STORE gvNm(:,:,k,bi,bj,2)
516    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
517    #  endif
518    #  ifdef ALLOW_CD_CODE
519    CADJ STORE uNM1(:,:,k,bi,bj)
520    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
521    CADJ STORE vNM1(:,:,k,bi,bj)
522    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
523    CADJ STORE uVelD(:,:,k,bi,bj)
524    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
525    CADJ STORE vVelD(:,:,k,bi,bj)
526    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
527    #  endif
528    # endif /* NONLIN_FRSURF */
529    #endif /* ALLOW_AUTODIFF_TAMC */
530    
531    C--      Integrate hydrostatic balance for phiHyd with BC of phiHyd(z=0)=0
532             IF ( implicitIntGravWave ) THEN
533             CALL CALC_PHI_HYD(             CALL CALC_PHI_HYD(
534       I        bi,bj,iMin,iMax,jMin,jMax,k,       I        bi,bj,iMin,iMax,jMin,jMax,k,
535       I        gTnm1, gSnm1,       I        gT, gS,
536       U        phiHyd,       U        phiHydF,
537       I        myThid )       O        phiHydC, dPhiHydX, dPhiHydY,
538         I        myTime, myIter, myThid )
539           ELSE           ELSE
540             CALL CALC_PHI_HYD(             CALL CALC_PHI_HYD(
541       I        bi,bj,iMin,iMax,jMin,jMax,k,       I        bi,bj,iMin,iMax,jMin,jMax,k,
542       I        theta, salt,       I        theta, salt,
543       U        phiHyd,       U        phiHydF,
544       I        myThid )       O        phiHydC, dPhiHydX, dPhiHydY,
545         I        myTime, myIter, myThid )
546             ENDIF
547    #ifdef ALLOW_DIAGNOSTICS
548             IF ( dPhiHydDiagIsOn ) THEN
549               tmpFac = -1. _d 0
550               CALL DIAGNOSTICS_SCALE_FILL( dPhiHydX, tmpFac, 1,
551         &                           'Um_dPHdx', k, 1, 2, bi, bj, myThid )
552               CALL DIAGNOSTICS_SCALE_FILL( dPhiHydY, tmpFac, 1,
553         &                           'Vm_dPHdy', k, 1, 2, bi, bj, myThid )
554           ENDIF           ENDIF
555    #endif /* ALLOW_DIAGNOSTICS */
556    
557  C--      Calculate accelerations in the momentum equations (gU, gV, ...)  C--      Calculate accelerations in the momentum equations (gU, gV, ...)
558  C        and step forward storing the result in gUnm1, gVnm1, etc...  C        and step forward storing the result in gU, gV, etc...
559           IF ( momStepping ) THEN           IF ( momStepping ) THEN
560             CALL CALC_MOM_RHS(  #ifdef ALLOW_AUTODIFF
561       I         bi,bj,iMin,iMax,jMin,jMax,k,kup,kDown,             DO j=1-OLy,sNy+OLy
562       I         phiHyd,KappaRU,KappaRV,              DO i=1-OLx,sNx+OLx
563       U         fVerU, fVerV,                guDissip(i,j)  = 0. _d 0
564       I         myTime, myThid)                gvDissip(i,j)  = 0. _d 0
565                ENDDO
566               ENDDO
567    #endif /* ALLOW_AUTODIFF */
568    #ifdef ALLOW_AUTODIFF_TAMC
569    # if (defined NONLIN_FRSURF) && (defined ALLOW_MOM_FLUXFORM)
570    #  ifndef DISABLE_RSTAR_CODE
571    CADJ STORE dWtransC(:,:,bi,bj)
572    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
573    CADJ STORE dWtransU(:,:,bi,bj)
574    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
575    CADJ STORE dWtransV(:,:,bi,bj)
576    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
577    #  endif
578    # endif /* NONLIN_FRSURF and ALLOW_MOM_FLUXFORM */
579    # if (defined NONLIN_FRSURF) || (defined ALLOW_DEPTH_CONTROL)
580    CADJ STORE fVerU(:,:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
581    CADJ STORE fVerV(:,:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
582    # endif
583    #endif /* ALLOW_AUTODIFF_TAMC */
584               IF (.NOT. vectorInvariantMomentum) THEN
585    #ifdef ALLOW_MOM_FLUXFORM
586                  CALL MOM_FLUXFORM(
587         I         bi,bj,k,iMin,iMax,jMin,jMax,
588         I         KappaRU, KappaRV,
589         U         fVerU(1-OLx,1-OLy,kUp),   fVerV(1-OLx,1-OLy,kUp),
590         O         fVerU(1-OLx,1-OLy,kDown), fVerV(1-OLx,1-OLy,kDown),
591         O         guDissip, gvDissip,
592         I         myTime, myIter, myThid)
593    #endif
594               ELSE
595    #ifdef ALLOW_MOM_VECINV
596                 CALL MOM_VECINV(
597         I         bi,bj,k,iMin,iMax,jMin,jMax,
598         I         KappaRU, KappaRV,
599         I         fVerU(1-OLx,1-OLy,kUp),   fVerV(1-OLx,1-OLy,kUp),
600         O         fVerU(1-OLx,1-OLy,kDown), fVerV(1-OLx,1-OLy,kDown),
601         O         guDissip, gvDissip,
602         I         myTime, myIter, myThid)
603    #endif
604               ENDIF
605    
606    #ifdef ALLOW_SMAG_3D
607               IF ( useSmag3D ) THEN
608                 CALL MOM_CALC_SMAG_3D(
609         I         str11, str22, str33, str12, str13, str23,
610         O         viscAh3d_00, viscAh3d_12, viscAh3d_13, viscAh3d_23,
611         I         smag3D_hLsC, smag3D_hLsW, smag3D_hLsS, smag3D_hLsZ,
612         I         k, bi, bj, myThid )
613                 CALL MOM_UV_SMAG_3D(
614         I         str11, str22, str12, str13, str23,
615         I         viscAh3d_00, viscAh3d_12, viscAh3d_13, viscAh3d_23,
616         O         addDissU, addDissV,
617         I         iMin,iMax,jMin,jMax, k, bi, bj, myThid )
618                 DO j= jMin,jMax
619                  DO i= iMin,iMax
620                   guDissip(i,j) = guDissip(i,j) + addDissU(i,j)
621                   gvDissip(i,j) = gvDissip(i,j) + addDissV(i,j)
622                  ENDDO
623                 ENDDO
624               ENDIF
625    #endif /* ALLOW_SMAG_3D */
626    
627             CALL TIMESTEP(             CALL TIMESTEP(
628       I         bi,bj,iMin,iMax,jMin,jMax,k,       I         bi,bj,iMin,iMax,jMin,jMax,k,
629       I         phiHyd, phiSurfX, phiSurfY,       I         dPhiHydX,dPhiHydY, phiSurfX, phiSurfY,
630       I         myIter, myThid)       I         guDissip, gvDissip,
631         I         myTime, myIter, myThid)
632    
 #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 */  
633           ENDIF           ENDIF
634    
   
635  C--     end of dynamics k loop (1:Nr)  C--     end of dynamics k loop (1:Nr)
636          ENDDO          ENDDO
637    
638    C--     Implicit Vertical advection & viscosity
639    #if (defined (INCLUDE_IMPLVERTADV_CODE) && \
640  C--     Implicit viscosity       defined (ALLOW_MOM_COMMON) && !(defined ALLOW_AUTODIFF_TAMC))
641          IF (implicitViscosity.AND.momStepping) THEN          IF ( momImplVertAdv ) THEN
642              CALL MOM_U_IMPLICIT_R( kappaRU,
643         I                           bi, bj, myTime, myIter, myThid )
644              CALL MOM_V_IMPLICIT_R( kappaRV,
645         I                           bi, bj, myTime, myIter, myThid )
646            ELSEIF ( implicitViscosity ) THEN
647    #else /* INCLUDE_IMPLVERTADV_CODE */
648            IF     ( implicitViscosity ) THEN
649    #endif /* INCLUDE_IMPLVERTADV_CODE */
650  #ifdef    ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
651            idkey = iikey + 3  CADJ STORE gU(:,:,:,bi,bj) = comlev1_bibj , key=idynkey, byte=isbyte
 CADJ STORE gUNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  
652  #endif    /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
653            CALL IMPLDIFF(            CALL IMPLDIFF(
654       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
655       I         deltaTmom, KappaRU,recip_HFacW,       I         -1, KappaRU, recip_hFacW(1-OLx,1-OLy,1,bi,bj),
656       U         gUNm1,       U         gU,
657       I         myThid )       I         myThid )
658  #ifdef    ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
659            idkey = iikey + 4  CADJ STORE gV(:,:,:,bi,bj) = comlev1_bibj , key=idynkey, byte=isbyte
 CADJ STORE gVNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  
660  #endif    /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
661            CALL IMPLDIFF(            CALL IMPLDIFF(
662       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
663       I         deltaTmom, KappaRV,recip_HFacS,       I         -2, KappaRV, recip_hFacS(1-OLx,1-OLy,1,bi,bj),
664       U         gVNm1,       U         gV,
665       I         myThid )       I         myThid )
666            ENDIF
667    
668  #ifdef   ALLOW_OBCS  #ifdef ALLOW_OBCS
669  C--      Apply open boundary conditions  C--      Apply open boundary conditions
670           IF (useOBCS) THEN          IF ( useOBCS ) THEN
671             DO K=1,Nr  C--      but first save intermediate velocities to be used in the
672               CALL OBCS_APPLY_UV( bi, bj, k, gUnm1, gVnm1, myThid )  C        next time step for the Stevens boundary conditions
673             ENDDO            CALL OBCS_SAVE_UV_N(
674           END IF       I        bi, bj, iMin, iMax, jMin, jMax, 0,
675  #endif   /* ALLOW_OBCS */       I        gU, gV, myThid )
676              CALL OBCS_APPLY_UV( bi, bj, 0, gU, gV, myThid )
677            ENDIF
678    #endif /* ALLOW_OBCS */
679    
680  #ifdef    INCLUDE_CD_CODE  #ifdef    ALLOW_CD_CODE
681            IF (implicitViscosity.AND.useCDscheme) THEN
682  #ifdef    ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
683            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  
684  #endif    /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
685            CALL IMPLDIFF(            CALL IMPLDIFF(
686       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
687       I         deltaTmom, KappaRU,recip_HFacW,       I         0, KappaRU, recip_hFacW(1-OLx,1-OLy,1,bi,bj),
688       U         vVelD,       U         vVelD,
689       I         myThid )       I         myThid )
690  #ifdef    ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
691            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  
692  #endif    /* ALLOW_AUTODIFF_TAMC */  #endif    /* ALLOW_AUTODIFF_TAMC */
693            CALL IMPLDIFF(            CALL IMPLDIFF(
694       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
695       I         deltaTmom, KappaRV,recip_HFacS,       I         0, KappaRV, recip_hFacS(1-OLx,1-OLy,1,bi,bj),
696       U         uVelD,       U         uVelD,
697       I         myThid )       I         myThid )
 #endif    /* INCLUDE_CD_CODE */  
 C--     End If implicitViscosity.AND.momStepping  
698          ENDIF          ENDIF
699    #endif    /* ALLOW_CD_CODE */
700  Cjmc : add for phiHyd output <- but not working if multi tile per CPU  C--     End implicit Vertical advection & viscosity
701  c       IF ( DIFFERENT_MULTIPLE(dumpFreq,myTime+deltaTClock,myTime)  
702  c    &  .AND. buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
703  c         WRITE(suff,'(I10.10)') myIter+1  
704  c         CALL WRITE_FLD_XYZ_RL('PH.',suff,phiHyd,myIter+1,myThid)  #ifdef ALLOW_NONHYDROSTATIC
705  c       ENDIF  C--   Step forward W field in N-H algorithm
706  Cjmc(end)          IF ( nonHydrostatic ) THEN
707    #ifdef ALLOW_DEBUG
708  #ifdef ALLOW_TIMEAVE           IF (debugMode) CALL DEBUG_CALL('CALC_GW', myThid )
709          IF (taveFreq.GT.0.) THEN  #endif
710            CALL TIMEAVE_CUMULATE(phiHydtave, phiHyd, Nr,           CALL TIMER_START('CALC_GW          [DYNAMICS]',myThid)
711       I                              deltaTclock, bi, bj, myThid)           CALL CALC_GW(
712            IF (ivdc_kappa.NE.0.) THEN       I                 bi,bj, KappaRU, KappaRV,
713              CALL TIMEAVE_CUMULATE(ConvectCountTave, ConvectCount, Nr,       I                 str13, str23, str33,
714       I                              deltaTclock, bi, bj, myThid)       I                 viscAh3d_00, viscAh3d_13, viscAh3d_23,
715            ENDIF       I                 myTime, myIter, myThid )
716          ENDIF          ENDIF
717  #endif /* ALLOW_TIMEAVE */          IF ( nonHydrostatic.OR.implicitIntGravWave )
718         &   CALL TIMESTEP_WVEL( bi,bj, myTime, myIter, myThid )
719            IF ( nonHydrostatic )
720         &   CALL TIMER_STOP ('CALC_GW          [DYNAMICS]',myThid)
721    #endif
722    
723    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
724    
725    C-    end of bi,bj loops
726         ENDDO         ENDDO
727        ENDDO        ENDDO
728    
729    #ifdef ALLOW_OBCS
730          IF (useOBCS) THEN
731            CALL OBCS_EXCHANGES( myThid )
732          ENDIF
733    #endif
734    
735    Cml(
736    C     In order to compare the variance of phiHydLow of a p/z-coordinate
737    C     run with etaH of a z/p-coordinate run the drift of phiHydLow
738    C     has to be removed by something like the following subroutine:
739    C      CALL REMOVE_MEAN_RL( 1, phiHydLow, maskInC, maskInC, rA, drF,
740    C     &                     'phiHydLow', myTime, myThid )
741    Cml)
742    
743    #ifdef ALLOW_DIAGNOSTICS
744          IF ( useDiagnostics ) THEN
745    
746           CALL DIAGNOSTICS_FILL(totPhihyd,'PHIHYD  ',0,Nr,0,1,1,myThid)
747           CALL DIAGNOSTICS_FILL(phiHydLow,'PHIBOT  ',0, 1,0,1,1,myThid)
748    
749           tmpFac = 1. _d 0
750           CALL DIAGNOSTICS_SCALE_FILL(totPhihyd,tmpFac,2,
751         &                                 'PHIHYDSQ',0,Nr,0,1,1,myThid)
752    
753           CALL DIAGNOSTICS_SCALE_FILL(phiHydLow,tmpFac,2,
754         &                                 'PHIBOTSQ',0, 1,0,1,1,myThid)
755    
756          ENDIF
757    #endif /* ALLOW_DIAGNOSTICS */
758    
759    #ifdef ALLOW_DEBUG
760          IF ( debugLevel .GE. debLevD ) THEN
761           CALL DEBUG_STATS_RL(1,EtaN,'EtaN (DYNAMICS)',myThid)
762           CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (DYNAMICS)',myThid)
763           CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (DYNAMICS)',myThid)
764           CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (DYNAMICS)',myThid)
765           CALL DEBUG_STATS_RL(Nr,theta,'Theta (DYNAMICS)',myThid)
766           CALL DEBUG_STATS_RL(Nr,salt,'Salt (DYNAMICS)',myThid)
767           CALL DEBUG_STATS_RL(Nr,gU,'Gu (DYNAMICS)',myThid)
768           CALL DEBUG_STATS_RL(Nr,gV,'Gv (DYNAMICS)',myThid)
769           CALL DEBUG_STATS_RL(Nr,gT,'Gt (DYNAMICS)',myThid)
770           CALL DEBUG_STATS_RL(Nr,gS,'Gs (DYNAMICS)',myThid)
771    #ifndef ALLOW_ADAMSBASHFORTH_3
772           CALL DEBUG_STATS_RL(Nr,guNm1,'GuNm1 (DYNAMICS)',myThid)
773           CALL DEBUG_STATS_RL(Nr,gvNm1,'GvNm1 (DYNAMICS)',myThid)
774           CALL DEBUG_STATS_RL(Nr,gtNm1,'GtNm1 (DYNAMICS)',myThid)
775           CALL DEBUG_STATS_RL(Nr,gsNm1,'GsNm1 (DYNAMICS)',myThid)
776    #endif
777          ENDIF
778    #endif
779    
780    #ifdef DYNAMICS_GUGV_EXCH_CHECK
781    C- jmc: For safety checking only: This Exchange here should not change
782    C       the solution. If solution changes, it means something is wrong,
783    C       but it does not mean that it is less wrong with this exchange.
784          IF ( debugLevel .GE. debLevE ) THEN
785           CALL EXCH_UV_XYZ_RL(gU,gV,.TRUE.,myThid)
786          ENDIF
787    #endif
788    
789    #ifdef ALLOW_DEBUG
790          IF (debugMode) CALL DEBUG_LEAVE( 'DYNAMICS', myThid )
791    #endif
792    
793        RETURN        RETURN
794        END        END
Legend:
Removed from v.1.67  
changed lines
  Added in v.1.168
  ViewVC Help
Powered by ViewVC 1.1.22