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