/[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.20 by cnh, Mon Jun 15 05:17:42 1998 UTC revision 1.74 by heimbach, Mon Jul 30 20:37:45 2001 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    C $Name$
3    
4  #include "CPP_EEOPTIONS.h"  #include "CPP_OPTIONS.h"
5    
6        SUBROUTINE DYNAMICS(myTime, myIter, myThid)        SUBROUTINE DYNAMICS(myTime, myIter, myThid)
7  C     /==========================================================\  C     /==========================================================\
# Line 20  C     | ===== Line 21  C     | =====
21  C     | C*P* comments indicating place holders for which code is |  C     | C*P* comments indicating place holders for which code is |
22  C     |      presently being developed.                          |  C     |      presently being developed.                          |
23  C     \==========================================================/  C     \==========================================================/
24          IMPLICIT NONE
25    
26  C     == Global variables ===  C     == Global variables ===
27  #include "SIZE.h"  #include "SIZE.h"
28  #include "EEPARAMS.h"  #include "EEPARAMS.h"
 #include "CG2D.h"  
29  #include "PARAMS.h"  #include "PARAMS.h"
30  #include "DYNVARS.h"  #include "DYNVARS.h"
31    #include "GRID.h"
32    #ifdef ALLOW_PASSIVE_TRACER
33    #include "TR1.h"
34    #endif
35    
36    #ifdef ALLOW_AUTODIFF_TAMC
37    # include "tamc.h"
38    # include "tamc_keys.h"
39    # include "FFIELDS.h"
40    # ifdef ALLOW_KPP
41    #  include "KPP.h"
42    # endif
43    # ifdef ALLOW_GMREDI
44    #  include "GMREDI.h"
45    # endif
46    #endif /* ALLOW_AUTODIFF_TAMC */
47    
48    #ifdef ALLOW_TIMEAVE
49    #include "TIMEAVE_STATV.h"
50    #endif
51    
52  C     == Routine arguments ==  C     == Routine arguments ==
53  C     myTime - Current time in simulation  C     myTime - Current time in simulation
54  C     myIter - Current iteration number in simulation  C     myIter - Current iteration number in simulation
55  C     myThid - Thread number for this instance of the routine.  C     myThid - Thread number for this instance of the routine.
       INTEGER myThid  
56        _RL myTime        _RL myTime
57        INTEGER myIter        INTEGER myIter
58          INTEGER myThid
59    
60  C     == Local variables  C     == Local variables
61  C     xA, yA                 - Per block temporaries holding face areas  C     xA, yA                 - Per block temporaries holding face areas
62  C     uTrans, vTrans, wTrans - Per block temporaries holding flow transport  C     uTrans, vTrans, rTrans - Per block temporaries holding flow
63  C     wVel                     o uTrans: Zonal transport  C                              transport
64    C                              o uTrans: Zonal transport
65  C                              o vTrans: Meridional transport  C                              o vTrans: Meridional transport
66  C                              o wTrans: Vertical transport  C                              o rTrans: Vertical transport
67  C                              o wVel:   Vertical velocity at upper and lower  C     maskUp                   o maskUp: land/water mask for W points
68  C                                        cell faces.  C     fVer[STUV]               o fVer: Vertical flux term - note fVer
 C     maskC,maskUp             o maskC: land/water mask for tracer cells  
 C                              o maskUp: land/water mask for W points  
 C     aTerm, xTerm, cTerm    - Work arrays for holding separate terms in  
 C     mTerm, pTerm,            tendency equations.  
 C     fZon, fMer, fVer[STUV]   o aTerm: Advection term  
 C                              o xTerm: Mixing term  
 C                              o cTerm: Coriolis term  
 C                              o mTerm: Metric term  
 C                              o pTerm: Pressure term  
 C                              o fZon: Zonal flux term  
 C                              o fMer: Meridional flux term  
 C                              o fVer: Vertical flux term - note fVer  
69  C                                      is "pipelined" in the vertical  C                                      is "pipelined" in the vertical
70  C                                      so we need an fVer for each  C                                      so we need an fVer for each
71  C                                      variable.  C                                      variable.
72  C     iMin, iMax - Ranges and sub-block indices on which calculations  C     rhoK, rhoKM1   - Density at current level, and level above
73  C     jMin, jMax   are applied.  C     phiHyd         - Hydrostatic part of the potential phiHydi.
74    C                      In z coords phiHydiHyd is the hydrostatic
75    C                      Potential (=pressure/rho0) anomaly
76    C                      In p coords phiHydiHyd is the geopotential
77    C                      surface height anomaly.
78    C     phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean)
79    C     phiSurfY             or geopotentiel (atmos) in X and Y direction
80    C     KappaRT,       - Total diffusion in vertical for T and S.
81    C     KappaRS          (background + spatially varying, isopycnal term).
82    C     iMin, iMax     - Ranges and sub-block indices on which calculations
83    C     jMin, jMax       are applied.
84  C     bi, bj  C     bi, bj
85  C     k, kUp, kDown, kM1 - Index for layer above and below. kUp and kDown  C     k, kup,        - Index for layer above and below. kup and kDown
86  C                          are switched with layer to be the appropriate index  C     kDown, km1       are switched with layer to be the appropriate
87  C                          into fVerTerm  C                      index into fVerTerm.
88        _RS xA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  C     tauAB - Adams-Bashforth timestepping weight: 0=forward ; 1/2=Adams-Bashf.
89        _RS yA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
90        _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
91        _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
92        _RL wTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
93        _RL wVel  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
94        _RS maskC (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
95        _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
96        _RL aTerm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fVerS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
97        _RL xTerm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fVerTr1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
98        _RL cTerm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fVerU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
99        _RL mTerm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fVerV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
100        _RL pTerm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiHyd  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
101        _RL fZon  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
102        _RL fMer  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
103        _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
104        _RL fVerS (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
105        _RL fVerU (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
106        _RL fVerV (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
107        _RL pH    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nz)        _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
108        _RL rhokm1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL KappaRV (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
109        _RL rhokp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
110        _RL rhok  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
111        _RL rhotmp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
112        _RL pSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL tauAB
113        _RL pSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
114        _RL K13   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nz)  C This is currently used by IVDC and Diagnostics
115        _RL K23   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nz)        _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
       _RL K33   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nz)  
       _RL KapGM (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL KappaZT(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nz)  
       _RL KappaZS(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nz)  
116    
117        INTEGER iMin, iMax        INTEGER iMin, iMax
118        INTEGER jMin, jMax        INTEGER jMin, jMax
119        INTEGER bi, bj        INTEGER bi, bj
120        INTEGER i, j        INTEGER i, j
121        INTEGER k, kM1, kUp, kDown        INTEGER k, km1, kup, kDown
       LOGICAL BOTTOM_LAYER  
122    
123    Cjmc : add for phiHyd output <- but not working if multi tile per CPU
124    c     CHARACTER*(MAX_LEN_MBUF) suff
125    c     LOGICAL  DIFFERENT_MULTIPLE
126    c     EXTERNAL DIFFERENT_MULTIPLE
127    Cjmc(end)
128    
129  C---    The algorithm...  C---    The algorithm...
130  C  C
131  C       "Correction Step"  C       "Correction Step"
# Line 119  C       "Calculation of Gs" Line 140  C       "Calculation of Gs"
140  C       ===================  C       ===================
141  C       This is where all the accelerations and tendencies (ie.  C       This is where all the accelerations and tendencies (ie.
142  C       physics, parameterizations etc...) are calculated  C       physics, parameterizations etc...) are calculated
 C         w = sum_z ( div. u[n] )  
143  C         rho = rho ( theta[n], salt[n] )  C         rho = rho ( theta[n], salt[n] )
144    C         b   = b(rho, theta)
145  C         K31 = K31 ( rho )  C         K31 = K31 ( rho )
146  C         Gu[n] = Gu( u[n], v[n], w, rho, Ph, ... )  C         Gu[n] = Gu( u[n], v[n], wVel, b, ... )
147  C         Gv[n] = Gv( u[n], v[n], w, rho, Ph, ... )  C         Gv[n] = Gv( u[n], v[n], wVel, b, ... )
148  C         Gt[n] = Gt( theta[n], u[n], v[n], w, K31, ... )  C         Gt[n] = Gt( theta[n], u[n], v[n], wVel, K31, ... )
149  C         Gs[n] = Gs( salt[n], u[n], v[n], w, K31, ... )  C         Gs[n] = Gs( salt[n], u[n], v[n], wVel, K31, ... )
150  C  C
151  C       "Time-stepping" or "Prediction"  C       "Time-stepping" or "Prediction"
152  C       ================================  C       ================================
# Line 149  C         (1 + dt * K * d_zz) theta[n] = Line 170  C         (1 + dt * K * d_zz) theta[n] =
170  C         (1 + dt * K * d_zz) salt[n] = salt*  C         (1 + dt * K * d_zz) salt[n] = salt*
171  C---  C---
172    
173    #ifdef ALLOW_AUTODIFF_TAMC
174    C--   dummy statement to end declaration part
175          ikey = 1
176    #endif /* ALLOW_AUTODIFF_TAMC */
177    
178  C--   Set up work arrays with valid (i.e. not NaN) values  C--   Set up work arrays with valid (i.e. not NaN) values
179  C     These inital values do not alter the numerical results. They  C     These inital values do not alter the numerical results. They
180  C     just ensure that all memory references are to valid floating  C     just ensure that all memory references are to valid floating
# Line 160  C     uninitialised but inert locations. Line 186  C     uninitialised but inert locations.
186          yA(i,j)      = 0. _d 0          yA(i,j)      = 0. _d 0
187          uTrans(i,j)  = 0. _d 0          uTrans(i,j)  = 0. _d 0
188          vTrans(i,j)  = 0. _d 0          vTrans(i,j)  = 0. _d 0
189          aTerm(i,j)   = 0. _d 0          DO k=1,Nr
190          xTerm(i,j)   = 0. _d 0           phiHyd(i,j,k)  = 0. _d 0
191          cTerm(i,j)   = 0. _d 0           KappaRU(i,j,k) = 0. _d 0
192          mTerm(i,j)   = 0. _d 0           KappaRV(i,j,k) = 0. _d 0
193          pTerm(i,j)   = 0. _d 0           sigmaX(i,j,k) = 0. _d 0
194          fZon(i,j)    = 0. _d 0           sigmaY(i,j,k) = 0. _d 0
195          fMer(i,j)    = 0. _d 0           sigmaR(i,j,k) = 0. _d 0
         DO K=1,nZ  
          pH (i,j,k)  = 0. _d 0  
          K13(i,j,k) = 0. _d 0  
          K23(i,j,k) = 0. _d 0  
          K33(i,j,k) = 0. _d 0  
          KappaZT(i,j,k) = 0. _d 0  
196          ENDDO          ENDDO
197          rhokm1(i,j)  = 0. _d 0          rhoKM1 (i,j) = 0. _d 0
198          rhok  (i,j)  = 0. _d 0          rhok   (i,j) = 0. _d 0
199          rhokp1(i,j)  = 0. _d 0          phiSurfX(i,j) = 0. _d 0
200          rhotmp(i,j)  = 0. _d 0          phiSurfY(i,j) = 0. _d 0
         maskC (i,j)  = 0. _d 0  
201         ENDDO         ENDDO
202        ENDDO        ENDDO
203    
204    
205    #ifdef ALLOW_AUTODIFF_TAMC
206    C--   HPF directive to help TAMC
207    CHPF$ INDEPENDENT
208    #endif /* ALLOW_AUTODIFF_TAMC */
209    
210        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
211    
212    #ifdef ALLOW_AUTODIFF_TAMC
213    C--    HPF directive to help TAMC
214    CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS,fVerU,fVerV
215    CHPF$&                  ,phiHyd,utrans,vtrans,xA,yA
216    CHPF$&                  ,KappaRT,KappaRS,KappaRU,KappaRV
217    CHPF$&                  )
218    #endif /* ALLOW_AUTODIFF_TAMC */
219    
220         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
221    
222    #ifdef ALLOW_AUTODIFF_TAMC
223              act1 = bi - myBxLo(myThid)
224              max1 = myBxHi(myThid) - myBxLo(myThid) + 1
225    
226              act2 = bj - myByLo(myThid)
227              max2 = myByHi(myThid) - myByLo(myThid) + 1
228    
229              act3 = myThid - 1
230              max3 = nTx*nTy
231    
232              act4 = ikey_dynamics - 1
233    
234              ikey = (act1 + 1) + act2*max1
235         &                      + act3*max1*max2
236         &                      + act4*max1*max2*max3
237    #endif /* ALLOW_AUTODIFF_TAMC */
238    
239  C--     Set up work arrays that need valid initial values  C--     Set up work arrays that need valid initial values
240          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
241           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
242            wTrans(i,j)  = 0. _d 0            rTrans (i,j)   = 0. _d 0
243            wVel  (i,j,1) = 0. _d 0            fVerT  (i,j,1) = 0. _d 0
244            wVel  (i,j,2) = 0. _d 0            fVerT  (i,j,2) = 0. _d 0
245            fVerT(i,j,1) = 0. _d 0            fVerS  (i,j,1) = 0. _d 0
246            fVerT(i,j,2) = 0. _d 0            fVerS  (i,j,2) = 0. _d 0
247            fVerS(i,j,1) = 0. _d 0            fVerTr1(i,j,1) = 0. _d 0
248            fVerS(i,j,2) = 0. _d 0            fVerTr1(i,j,2) = 0. _d 0
249            fVerU(i,j,1) = 0. _d 0            fVerU  (i,j,1) = 0. _d 0
250            fVerU(i,j,2) = 0. _d 0            fVerU  (i,j,2) = 0. _d 0
251            fVerV(i,j,1) = 0. _d 0            fVerV  (i,j,1) = 0. _d 0
252            fVerV(i,j,2) = 0. _d 0            fVerV  (i,j,2) = 0. _d 0
253            pH(i,j,1) = 0. _d 0           ENDDO
254            K13(i,j,1) = 0. _d 0          ENDDO
255            K23(i,j,1) = 0. _d 0  
256            K33(i,j,1) = 0. _d 0          DO k=1,Nr
257            KapGM(i,j) = 0. _d 0           DO j=1-OLy,sNy+OLy
258              DO i=1-OLx,sNx+OLx
259    C This is currently also used by IVDC and Diagnostics
260               ConvectCount(i,j,k) = 0.
261               KappaRT(i,j,k) = 0. _d 0
262               KappaRS(i,j,k) = 0. _d 0
263              ENDDO
264           ENDDO           ENDDO
265          ENDDO          ENDDO
266    
# Line 212  C--     Set up work arrays that need val Line 269  C--     Set up work arrays that need val
269          jMin = 1-OLy+1          jMin = 1-OLy+1
270          jMax = sNy+OLy          jMax = sNy+OLy
271    
         K = 1  
         BOTTOM_LAYER = K .EQ. Nz  
272    
273  C--     Calculate gradient of surface pressure  #ifdef ALLOW_AUTODIFF_TAMC
274          CALL GRAD_PSURF(  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
275       I       bi,bj,iMin,iMax,jMin,jMax,  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
276       O       pSurfX,pSurfY,  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
277       I       myThid)  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
278    #ifdef ALLOW_PASSIVE_TRACER
279  C--     Update fields in top level according to tendency terms  CADJ STORE tr1  (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
280          CALL CORRECTION_STEP(  #endif
281       I       bi,bj,iMin,iMax,jMin,jMax,K,pSurfX,pSurfY,myThid)  #endif /* ALLOW_AUTODIFF_TAMC */
282    
283  C--     Density of 1st level (below W(1)) reference to level 1  C--     Start of diagnostic loop
284          CALL FIND_RHO(          DO k=Nr,1,-1
285       I     bi, bj, iMin, iMax, jMin, jMax, K, K, eosType,  
286       O     rhoKm1,  #ifdef ALLOW_AUTODIFF_TAMC
287       I     myThid )  C? Patrick, is this formula correct now that we change the loop range?
288    C? Do we still need this?
289          IF ( .NOT. BOTTOM_LAYER ) THEN  cph kkey formula corrected.
290  C--      Check static stability with layer below  cph Needed for rhok, rhokm1, in the case useGMREDI.
291  C        and mix as needed.           kkey = (ikey-1)*Nr + k
292           CALL FIND_RHO(  CADJ STORE rhokm1(:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte
293       I      bi, bj, iMin, iMax, jMin, jMax, K+1, K, eosType,  CADJ STORE rhok  (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte
294       O      rhoKp1,  #endif /* ALLOW_AUTODIFF_TAMC */
295       I      myThid )  
296           CALL CONVECT(  C--       Integrate continuity vertically for vertical velocity
297       I       bi,bj,iMin,iMax,jMin,jMax,K+1,rhoKm1,rhoKp1,            CALL INTEGRATE_FOR_W(
298       I       myTime,myIter,myThid)       I                         bi, bj, k, uVel, vVel,
299  C--      Recompute density after mixing       O                         wVel,
300           CALL FIND_RHO(       I                         myThid )
301       I      bi, bj, iMin, iMax, jMin, jMax, K, K, eosType,  
302       O      rhoKm1,  #ifdef    ALLOW_OBCS
303       I      myThid )  #ifdef    ALLOW_NONHYDROSTATIC
304          ENDIF  C--       Apply OBC to W if in N-H mode
305              IF (useOBCS.AND.nonHydrostatic) THEN
306                CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )
307              ENDIF
308    #endif    /* ALLOW_NONHYDROSTATIC */
309    #endif    /* ALLOW_OBCS */
310    
311    C--       Calculate gradients of potential density for isoneutral
312    C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)
313    c         IF ( k.GT.1 .AND. (useGMRedi.OR.ivdc_kappa.NE.0.) ) THEN
314              IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.) ) THEN
315    #ifdef ALLOW_AUTODIFF_TAMC
316    CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
317    CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
318    #endif /* ALLOW_AUTODIFF_TAMC */
319                CALL FIND_RHO(
320         I        bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,
321         I        theta, salt,
322         O        rhoK,
323         I        myThid )
324                IF (k.GT.1) THEN
325    #ifdef ALLOW_AUTODIFF_TAMC
326    CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
327    CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
328    #endif /* ALLOW_AUTODIFF_TAMC */
329                 CALL FIND_RHO(
330         I        bi, bj, iMin, iMax, jMin, jMax, k-1, k, eosType,
331         I        theta, salt,
332         O        rhoKm1,
333         I        myThid )
334                ENDIF
335                CALL GRAD_SIGMA(
336         I             bi, bj, iMin, iMax, jMin, jMax, k,
337         I             rhoK, rhoKm1, rhoK,
338         O             sigmaX, sigmaY, sigmaR,
339         I             myThid )
340              ENDIF
341    
342  C--     Integrate hydrostatic balance for pH with BC of pH(z=0)=0  C--       Implicit Vertical Diffusion for Convection
343          CALL CALC_PH(  c ==> should use sigmaR !!!
344       I      bi,bj,iMin,iMax,jMin,jMax,K,rhoKm1,rhoKm1,            IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN
345       U      pH,              CALL CALC_IVDC(
346       I      myThid )       I        bi, bj, iMin, iMax, jMin, jMax, k,
347         I        rhoKm1, rhoK,
348         U        ConvectCount, KappaRT, KappaRS,
349         I        myTime, myIter, myThid)
350              ENDIF
351    
352          DO K=2,Nz  C--     end of diagnostic k loop (Nr:1)
353            ENDDO
354    
355           BOTTOM_LAYER = K .EQ. Nz  #ifdef ALLOW_AUTODIFF_TAMC
356    cph avoids recomputation of integrate_for_w
357    CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
358    #endif /* ALLOW_AUTODIFF_TAMC */
359    
360    #ifdef  ALLOW_OBCS
361    C--     Calculate future values on open boundaries
362            IF (useOBCS) THEN
363              CALL OBCS_CALC( bi, bj, myTime+deltaT,
364         I            uVel, vVel, wVel, theta, salt,
365         I            myThid )
366            ENDIF
367    #endif  /* ALLOW_OBCS */
368    
369  C--      Update fields in Kth level according to tendency terms  C--     Determines forcing terms based on external fields
370           CALL CORRECTION_STEP(  C       relaxation terms, etc.
371       I        bi,bj,iMin,iMax,jMin,jMax,K,pSurfX,pSurfY,myThid)          CALL EXTERNAL_FORCING_SURF(
372  C--      Density of K level (below W(K)) reference to K level       I             bi, bj, iMin, iMax, jMin, jMax,
373           CALL FIND_RHO(       I             myThid )
374       I      bi, bj, iMin, iMax, jMin, jMax,  K, K, eosType,  #ifdef ALLOW_AUTODIFF_TAMC
375       O      rhoK,  cph needed for KPP
376       I      myThid )  CADJ STORE surfacetendencyU(:,:,bi,bj)
377           IF ( .NOT. BOTTOM_LAYER ) THEN  CADJ &     = comlev1_bibj, key=ikey, byte=isbyte
378  C--       Check static stability with layer below  CADJ STORE surfacetendencyV(:,:,bi,bj)
379  C         and mix as needed.  CADJ &     = comlev1_bibj, key=ikey, byte=isbyte
380  C--       Density of K+1 level (below W(K+1)) reference to K level  CADJ STORE surfacetendencyS(:,:,bi,bj)
381            CALL FIND_RHO(  CADJ &     = comlev1_bibj, key=ikey, byte=isbyte
382       I       bi, bj, iMin, iMax, jMin, jMax,  K+1, K, eosType,  CADJ STORE surfacetendencyT(:,:,bi,bj)
383       O       rhoKp1,  CADJ &     = comlev1_bibj, key=ikey, byte=isbyte
384       I       myThid )  #endif /* ALLOW_AUTODIFF_TAMC */
385            CALL CONVECT(  
386       I        bi,bj,iMin,iMax,jMin,jMax,K+1,rhoK,rhoKp1,  #ifdef  ALLOW_GMREDI
387       I        myTime,myIter,myThid)  
388  C--       Recompute density after mixing  #ifdef ALLOW_AUTODIFF_TAMC
389            CALL FIND_RHO(  CADJ STORE sigmaX(:,:,:) = comlev1, key=ikey, byte=isbyte
390       I       bi, bj, iMin, iMax, jMin, jMax, K, K, eosType,  CADJ STORE sigmaY(:,:,:) = comlev1, key=ikey, byte=isbyte
391       O       rhoK,  CADJ STORE sigmaR(:,:,:) = comlev1, key=ikey, byte=isbyte
392       I       myThid )  #endif /* ALLOW_AUTODIFF_TAMC */
393           ENDIF  C--     Calculate iso-neutral slopes for the GM/Redi parameterisation
394  C--      Integrate hydrostatic balance for pH with BC of pH(z=0)=0          IF (useGMRedi) THEN
395           CALL CALC_PH(            DO k=1,Nr
396       I       bi,bj,iMin,iMax,jMin,jMax,K,rhoKm1,rhoK,              CALL GMREDI_CALC_TENSOR(
397       U       pH,       I             bi, bj, iMin, iMax, jMin, jMax, k,
398       I       myThid )       I             sigmaX, sigmaY, sigmaR,
 C--      Calculate iso-neutral slopes for the GM/Redi parameterisation  
          CALL FIND_RHO(  
      I      bi, bj, iMin, iMax, jMin, jMax, K-1, K, eosType,  
      O      rhoTmp,  
      I      myThid )  
          CALL CALC_ISOSLOPES(  
      I             bi, bj, iMin, iMax, jMin, jMax, K,  
      I             rhoKm1, rhoK, rhotmp,  
      O             K13, K23, K33, KapGM,  
399       I             myThid )       I             myThid )
          DO J=jMin,jMax  
           DO I=iMin,iMax  
            rhoKm1(I,J)=rhoK(I,J)  
400            ENDDO            ENDDO
401           ENDDO  #ifdef ALLOW_AUTODIFF_TAMC
402            ELSE
403              DO k=1, Nr
404                CALL GMREDI_CALC_TENSOR_DUMMY(
405         I             bi, bj, iMin, iMax, jMin, jMax, k,
406         I             sigmaX, sigmaY, sigmaR,
407         I             myThid )
408              ENDDO
409    #endif /* ALLOW_AUTODIFF_TAMC */
410            ENDIF
411    
412          ENDDO ! K  #ifdef ALLOW_AUTODIFF_TAMC
413    CADJ STORE Kwx(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte
414    CADJ STORE Kwy(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte
415    CADJ STORE Kwz(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte
416    #endif /* ALLOW_AUTODIFF_TAMC */
417    
418    #endif  /* ALLOW_GMREDI */
419    
420    #ifdef  ALLOW_KPP
421    C--     Compute KPP mixing coefficients
422            IF (useKPP) THEN
423              CALL KPP_CALC(
424         I                  bi, bj, myTime, myThid )
425    #ifdef ALLOW_AUTODIFF_TAMC
426            ELSE
427              CALL KPP_CALC_DUMMY(
428         I                  bi, bj, myTime, myThid )
429    #endif /* ALLOW_AUTODIFF_TAMC */
430            ENDIF
431    
432          DO K = Nz, 1, -1  #ifdef ALLOW_AUTODIFF_TAMC
433           kM1  =max(1,k-1)   ! Points to level above k (=k-1)  CADJ STORE KPPghat   (:,:,:,bi,bj)
434           kUp  =1+MOD(k+1,2) ! Cycles through 1,2 to point to layer above  CADJ &   , KPPviscAz (:,:,:,bi,bj)
435           kDown=1+MOD(k,2)   ! Cycles through 2,1 to point to current layer  CADJ &   , KPPdiffKzT(:,:,:,bi,bj)
436           iMin = 1-OLx+2  CADJ &   , KPPdiffKzS(:,:,:,bi,bj)
437           iMax = sNx+OLx-1  CADJ &   , KPPfrac   (:,:  ,bi,bj)
438           jMin = 1-OLy+2  CADJ &                 = comlev1_bibj, key=ikey, byte=isbyte
439           jMax = sNy+OLy-1  #endif /* ALLOW_AUTODIFF_TAMC */
440    
441  C--      Get temporary terms used by tendency routines  #endif  /* ALLOW_KPP */
442           CALL CALC_COMMON_FACTORS (  
443       I        bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,  #ifdef ALLOW_AUTODIFF_TAMC
444       O        xA,yA,uTrans,vTrans,wTrans,wVel,maskC,maskUp,  CADJ STORE KappaRT(:,:,:)     = comlev1_bibj, key = ikey, byte = isbyte
445       I        myThid)  CADJ STORE KappaRS(:,:,:)     = comlev1_bibj, key = ikey, byte = isbyte
446    CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
447    CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
448    CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
449    CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
450    #ifdef ALLOW_PASSIVE_TRACER
451    CADJ STORE tr1  (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
452    #endif
453    #endif /* ALLOW_AUTODIFF_TAMC */
454    
455    #ifdef ALLOW_AIM
456    C       AIM - atmospheric intermediate model, physics package code.
457    C note(jmc) : phiHyd=0 at this point but is not really used in Molteni Physics
458            IF ( useAIM ) THEN
459             CALL TIMER_START('AIM_DO_ATMOS_PHYS      [DYNAMICS]', myThid)
460             CALL AIM_DO_ATMOS_PHYSICS( phiHyd, bi, bj, myTime, myThid )
461             CALL TIMER_STOP ('AIM_DO_ATMOS_PHYS      [DYNAMICS]', myThid)
462            ENDIF
463    #endif /* ALLOW_AIM */
464    
465    
466    C--     Start of thermodynamics loop
467            DO k=Nr,1,-1
468    #ifdef ALLOW_AUTODIFF_TAMC
469    C? Patrick Is this formula correct?
470    cph Yes, but I rewrote it.
471    cph Also, the KappaR? need the index and subscript k!
472             kkey = (ikey-1)*Nr + k
473    #endif /* ALLOW_AUTODIFF_TAMC */
474    
475    C--       km1    Points to level above k (=k-1)
476    C--       kup    Cycles through 1,2 to point to layer above
477    C--       kDown  Cycles through 2,1 to point to current layer
478    
479              km1  = MAX(1,k-1)
480              kup  = 1+MOD(k+1,2)
481              kDown= 1+MOD(k,2)
482    
483              iMin = 1-OLx
484              iMax = sNx+OLx
485              jMin = 1-OLy
486              jMax = sNy+OLy
487    
488    C--       Get temporary terms used by tendency routines
489              CALL CALC_COMMON_FACTORS (
490         I         bi,bj,iMin,iMax,jMin,jMax,k,
491         O         xA,yA,uTrans,vTrans,rTrans,maskUp,
492         I         myThid)
493    
494    #ifdef ALLOW_AUTODIFF_TAMC
495    CADJ STORE KappaRT(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte
496    CADJ STORE KappaRS(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte
497    #endif /* ALLOW_AUTODIFF_TAMC */
498    
499    #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL
500  C--      Calculate the total vertical diffusivity  C--      Calculate the total vertical diffusivity
501           CALL CALC_DIFFUSIVITY(           CALL CALC_DIFFUSIVITY(
502       I        bi,bj,iMin,iMax,jMin,jMax,K,       I        bi,bj,iMin,iMax,jMin,jMax,k,
503       I        maskC,maskUp,KapGM,K33,       I        maskUp,
504       O        KappaZT,KappaZS,       O        KappaRT,KappaRS,KappaRU,KappaRV,
505       I        myThid)       I        myThid)
506    #endif
507    
508  C--      Calculate accelerations in the momentum equations            iMin = 1-OLx+2
509           IF ( momStepping ) THEN            iMax = sNx+OLx-1
510            CALL CALC_MOM_RHS(            jMin = 1-OLy+2
511       I         bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,            jMax = sNy+OLy-1
      I         xA,yA,uTrans,vTrans,wTrans,wVel,maskC,  
      I         pH,  
      U         aTerm,xTerm,cTerm,mTerm,pTerm,  
      U         fZon, fMer, fVerU, fVerV,  
      I         myThid)  
          ENDIF  
512    
513  C--      Calculate active tracer tendencies  C--      Calculate active tracer tendencies (gT,gS,...)
514    C        and step forward storing result in gTnm1, gSnm1, etc.
515           IF ( tempStepping ) THEN           IF ( tempStepping ) THEN
516            CALL CALC_GT(             CALL CALC_GT(
517       I         bi,bj,iMin,iMax,jMin,jMax, k,kM1,kUp,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
518       I         xA,yA,uTrans,vTrans,wTrans,maskUp,       I         xA,yA,uTrans,vTrans,rTrans,maskUp,
519       I         K13,K23,KappaZT,KapGM,       I         KappaRT,
520       U         aTerm,xTerm,fZon,fMer,fVerT,       U         fVerT,
521       I         myThid)       I         myTime, myThid)
522               tauAB = 0.5d0 + abEps
523               CALL TIMESTEP_TRACER(
524         I         bi,bj,iMin,iMax,jMin,jMax,k,tauAB,
525         I         theta, gT,
526         U         gTnm1,
527         I         myIter, myThid)
528           ENDIF           ENDIF
529           IF ( saltStepping ) THEN           IF ( saltStepping ) THEN
530            CALL CALC_GS(             CALL CALC_GS(
531       I         bi,bj,iMin,iMax,jMin,jMax, k,kM1,kUp,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
532       I         xA,yA,uTrans,vTrans,wTrans,maskUp,       I         xA,yA,uTrans,vTrans,rTrans,maskUp,
533       I         K13,K23,KappaZS,KapGM,       I         KappaRS,
534       U         aTerm,xTerm,fZon,fMer,fVerS,       U         fVerS,
535       I         myThid)       I         myTime, myThid)
536               tauAB = 0.5d0 + abEps
537               CALL TIMESTEP_TRACER(
538         I         bi,bj,iMin,iMax,jMin,jMax,k,tauAB,
539         I         salt, gS,
540         U         gSnm1,
541         I         myIter, myThid)
542           ENDIF           ENDIF
543    #ifdef ALLOW_PASSIVE_TRACER
544             IF ( tr1Stepping ) THEN
545               CALL CALC_GTR1(
546         I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
547         I         xA,yA,uTrans,vTrans,rTrans,maskUp,
548         I         KappaRT,
549         U         fVerTr1,
550         I         myTime, myThid)
551               tauAB = 0.5d0 + abEps
552               CALL TIMESTEP_TRACER(
553         I         bi,bj,iMin,iMax,jMin,jMax,k,tauAB,
554         I         Tr1, gTr1,
555         U         gTr1NM1,
556         I         myIter, myThid)
557             ENDIF
558    #endif
559    
560    #ifdef   ALLOW_OBCS
561    C--      Apply open boundary conditions
562             IF (useOBCS) THEN
563               CALL OBCS_APPLY_TS( bi, bj, k, gTnm1, gSnm1, myThid )
564             END IF
565    #endif   /* ALLOW_OBCS */
566    
567    C--      Freeze water
568             IF (allowFreezing) THEN
569    #ifdef ALLOW_AUTODIFF_TAMC
570    CADJ STORE gTNm1(:,:,k,bi,bj) = comlev1_bibj_k
571    CADJ &   , key = kkey, byte = isbyte
572    #endif /* ALLOW_AUTODIFF_TAMC */
573                CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )
574             END IF
575    
576    C--     end of thermodynamic k loop (Nr:1)
577            ENDDO
578    
 C--      Prediction step (step forward all model variables)  
          CALL TIMESTEP(  
      I       bi,bj,iMin,iMax,jMin,jMax,K,  
      I       myThid)  
   
 C--      Diagnose barotropic divergence of predicted fields  
          CALL DIV_G(  
      I       bi,bj,iMin,iMax,jMin,jMax,K,  
      I       xA,yA,  
      I       myThid)  
579    
580          ENDDO ! K  #ifdef ALLOW_AUTODIFF_TAMC
581    C? Patrick? What about this one?
582    cph Keys iikey and idkey don't seem to be needed
583    cph since storing occurs on different tape for each
584    cph impldiff call anyways.
585    cph Thus, common block comlev1_impl isn't needed either.
586    cph Storing below needed in the case useGMREDI.
587            iikey = (ikey-1)*maximpl
588    #endif /* ALLOW_AUTODIFF_TAMC */
589    
590  C--     Implicit diffusion  C--     Implicit diffusion
591          IF (implicitDiffusion) THEN          IF (implicitDiffusion) THEN
592           CALL IMPLDIFF( bi, bj, iMin, iMax, jMin, jMax,  
593       I                  KappaZT,KappaZS,           IF (tempStepping) THEN
594       I                  myThid )  #ifdef ALLOW_AUTODIFF_TAMC
595                idkey = iikey + 1
596    CADJ STORE gTNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
597    #endif /* ALLOW_AUTODIFF_TAMC */
598                CALL IMPLDIFF(
599         I         bi, bj, iMin, iMax, jMin, jMax,
600         I         deltaTtracer, KappaRT, recip_HFacC,
601         U         gTNm1,
602         I         myThid )
603             ENDIF
604    
605             IF (saltStepping) THEN
606    #ifdef ALLOW_AUTODIFF_TAMC
607             idkey = iikey + 2
608    CADJ STORE gSNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
609    #endif /* ALLOW_AUTODIFF_TAMC */
610                CALL IMPLDIFF(
611         I         bi, bj, iMin, iMax, jMin, jMax,
612         I         deltaTtracer, KappaRS, recip_HFacC,
613         U         gSNm1,
614         I         myThid )
615             ENDIF
616    
617    #ifdef ALLOW_PASSIVE_TRACER
618             IF (tr1Stepping) THEN
619    #ifdef ALLOW_AUTODIFF_TAMC
620    CADJ STORE gTr1Nm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
621    #endif /* ALLOW_AUTODIFF_TAMC */
622              CALL IMPLDIFF(
623         I      bi, bj, iMin, iMax, jMin, jMax,
624         I      deltaTtracer, KappaRT, recip_HFacC,
625         U      gTr1Nm1,
626         I      myThid )
627             ENDIF
628    #endif
629    
630    #ifdef   ALLOW_OBCS
631    C--      Apply open boundary conditions
632             IF (useOBCS) THEN
633               DO K=1,Nr
634                 CALL OBCS_APPLY_TS( bi, bj, k, gTnm1, gSnm1, myThid )
635               ENDDO
636             END IF
637    #endif   /* ALLOW_OBCS */
638    
639    C--     End If implicitDiffusion
640            ENDIF
641    
642    C--     Start computation of dynamics
643            iMin = 1-OLx+2
644            iMax = sNx+OLx-1
645            jMin = 1-OLy+2
646            jMax = sNy+OLy-1
647    
648    C--     Explicit part of the Surface Potentiel Gradient (add in TIMESTEP)
649    C       (note: this loop will be replaced by CALL CALC_GRAD_ETA)
650            IF (implicSurfPress.NE.1.) THEN
651              CALL CALC_GRAD_PHI_SURF(
652         I         bi,bj,iMin,iMax,jMin,jMax,
653         I         etaN,
654         O         phiSurfX,phiSurfY,
655         I         myThid )                        
656            ENDIF
657    
658    C--     Start of dynamics loop
659            DO k=1,Nr
660    
661    C--       km1    Points to level above k (=k-1)
662    C--       kup    Cycles through 1,2 to point to layer above
663    C--       kDown  Cycles through 2,1 to point to current layer
664    
665              km1  = MAX(1,k-1)
666              kup  = 1+MOD(k+1,2)
667              kDown= 1+MOD(k,2)
668    
669    C--      Integrate hydrostatic balance for phiHyd with BC of
670    C        phiHyd(z=0)=0
671    C        distinguishe between Stagger and Non Stagger time stepping
672             IF (staggerTimeStep) THEN
673               CALL CALC_PHI_HYD(
674         I        bi,bj,iMin,iMax,jMin,jMax,k,
675         I        gTnm1, gSnm1,
676         U        phiHyd,
677         I        myThid )
678             ELSE
679               CALL CALC_PHI_HYD(
680         I        bi,bj,iMin,iMax,jMin,jMax,k,
681         I        theta, salt,
682         U        phiHyd,
683         I        myThid )
684             ENDIF
685    
686    C--      Calculate accelerations in the momentum equations (gU, gV, ...)
687    C        and step forward storing the result in gUnm1, gVnm1, etc...
688             IF ( momStepping ) THEN
689               CALL CALC_MOM_RHS(
690         I         bi,bj,iMin,iMax,jMin,jMax,k,kup,kDown,
691         I         phiHyd,KappaRU,KappaRV,
692         U         fVerU, fVerV,
693         I         myTime, myThid)
694               CALL TIMESTEP(
695         I         bi,bj,iMin,iMax,jMin,jMax,k,
696         I         phiHyd, phiSurfX, phiSurfY,
697         I         myIter, myThid)
698    
699    #ifdef   ALLOW_OBCS
700    C--      Apply open boundary conditions
701             IF (useOBCS) THEN
702               CALL OBCS_APPLY_UV( bi, bj, k, gUnm1, gVnm1, myThid )
703             END IF
704    #endif   /* ALLOW_OBCS */
705    
706    #ifdef   ALLOW_AUTODIFF_TAMC
707    #ifdef   INCLUDE_CD_CODE
708             ELSE
709               DO j=1-OLy,sNy+OLy
710                 DO i=1-OLx,sNx+OLx
711                   guCD(i,j,k,bi,bj) = 0.0
712                   gvCD(i,j,k,bi,bj) = 0.0
713                 END DO
714               END DO
715    #endif   /* INCLUDE_CD_CODE */
716    #endif   /* ALLOW_AUTODIFF_TAMC */
717             ENDIF
718    
719    
720    C--     end of dynamics k loop (1:Nr)
721            ENDDO
722    
723    
724    
725    C--     Implicit viscosity
726            IF (implicitViscosity.AND.momStepping) THEN
727    #ifdef    ALLOW_AUTODIFF_TAMC
728              idkey = iikey + 3
729    CADJ STORE gUNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
730    #endif    /* ALLOW_AUTODIFF_TAMC */
731              CALL IMPLDIFF(
732         I         bi, bj, iMin, iMax, jMin, jMax,
733         I         deltaTmom, KappaRU,recip_HFacW,
734         U         gUNm1,
735         I         myThid )
736    #ifdef    ALLOW_AUTODIFF_TAMC
737              idkey = iikey + 4
738    CADJ STORE gVNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
739    #endif    /* ALLOW_AUTODIFF_TAMC */
740              CALL IMPLDIFF(
741         I         bi, bj, iMin, iMax, jMin, jMax,
742         I         deltaTmom, KappaRV,recip_HFacS,
743         U         gVNm1,
744         I         myThid )
745    
746    #ifdef   ALLOW_OBCS
747    C--      Apply open boundary conditions
748             IF (useOBCS) THEN
749               DO K=1,Nr
750                 CALL OBCS_APPLY_UV( bi, bj, k, gUnm1, gVnm1, myThid )
751               ENDDO
752             END IF
753    #endif   /* ALLOW_OBCS */
754    
755    #ifdef    INCLUDE_CD_CODE
756    #ifdef    ALLOW_AUTODIFF_TAMC
757              idkey = iikey + 5
758    CADJ STORE vVelD(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
759    #endif    /* ALLOW_AUTODIFF_TAMC */
760              CALL IMPLDIFF(
761         I         bi, bj, iMin, iMax, jMin, jMax,
762         I         deltaTmom, KappaRU,recip_HFacW,
763         U         vVelD,
764         I         myThid )
765    #ifdef    ALLOW_AUTODIFF_TAMC
766              idkey = iikey + 6
767    CADJ STORE uVelD(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
768    #endif    /* ALLOW_AUTODIFF_TAMC */
769              CALL IMPLDIFF(
770         I         bi, bj, iMin, iMax, jMin, jMax,
771         I         deltaTmom, KappaRV,recip_HFacS,
772         U         uVelD,
773         I         myThid )
774    #endif    /* INCLUDE_CD_CODE */
775    C--     End If implicitViscosity.AND.momStepping
776          ENDIF          ENDIF
777    
778    Cjmc : add for phiHyd output <- but not working if multi tile per CPU
779    c       IF ( DIFFERENT_MULTIPLE(dumpFreq,myTime+deltaTClock,myTime)
780    c    &  .AND. buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN
781    c         WRITE(suff,'(I10.10)') myIter+1
782    c         CALL WRITE_FLD_XYZ_RL('PH.',suff,phiHyd,myIter+1,myThid)
783    c       ENDIF
784    Cjmc(end)
785    
786    #ifdef ALLOW_TIMEAVE
787            IF (taveFreq.GT.0.) THEN
788              CALL TIMEAVE_CUMUL_1T(phiHydtave, phiHyd, Nr,
789         I                              deltaTclock, bi, bj, myThid)
790              IF (ivdc_kappa.NE.0.) THEN
791                CALL TIMEAVE_CUMULATE(ConvectCountTave, ConvectCount, Nr,
792         I                              deltaTclock, bi, bj, myThid)
793              ENDIF
794            ENDIF
795    #endif /* ALLOW_TIMEAVE */
796    
797         ENDDO         ENDDO
798        ENDDO        ENDDO
799    
800  C     write(0,*) 'dynamics: pS ',minval(cg2d_x(1:sNx,1:sNy,:,:)),  #ifndef EXCLUDE_DEBUGMODE
801  C    &                           maxval(cg2d_x(1:sNx,1:sNy,:,:))        If (debugMode) THEN
802  C     write(0,*) 'dynamics: U  ',minval(uVel(1:sNx,1:sNy,1,:,:),mask=uVel(1:sNx,1:sNy,1,:,:).NE.0.),         CALL DEBUG_STATS_RL(1,EtaN,'EtaN (DYNAMICS)',myThid)
803  C    &                           maxval(uVel(1:sNx,1:sNy,1,:,:))         CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (DYNAMICS)',myThid)
804  C     write(0,*) 'dynamics: V  ',minval(vVel(1:sNx,1:sNy,1,:,:),mask=vVel(1:sNx,1:sNy,1,:,:).NE.0.),         CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (DYNAMICS)',myThid)
805  C    &                           maxval(vVel(1:sNx,1:sNy,1,:,:))         CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (DYNAMICS)',myThid)
806  C     write(0,*) 'dynamics: wVel(1) ',         CALL DEBUG_STATS_RL(Nr,theta,'Theta (DYNAMICS)',myThid)
807  C    &            minval(wVel(1:sNx,1:sNy,1),mask=wVel(1:sNx,1:sNy,1).NE.0.),         CALL DEBUG_STATS_RL(Nr,salt,'Salt (DYNAMICS)',myThid)
808  C    &            maxval(wVel(1:sNx,1:sNy,1))         CALL DEBUG_STATS_RL(Nr,Gu,'Gu (DYNAMICS)',myThid)
809  C     write(0,*) 'dynamics: wVel(2) ',         CALL DEBUG_STATS_RL(Nr,Gv,'Gv (DYNAMICS)',myThid)
810  C    &            minval(wVel(1:sNx,1:sNy,2),mask=wVel(1:sNx,1:sNy,2).NE.0.),         CALL DEBUG_STATS_RL(Nr,Gt,'Gt (DYNAMICS)',myThid)
811  C    &            maxval(wVel(1:sNx,1:sNy,2))         CALL DEBUG_STATS_RL(Nr,Gs,'Gs (DYNAMICS)',myThid)
812  cblk  write(0,*) 'dynamics: K13',minval(K13(1:sNx,1:sNy,:)),         CALL DEBUG_STATS_RL(Nr,GuNm1,'GuNm1 (DYNAMICS)',myThid)
813  cblk &                           maxval(K13(1:sNx,1:sNy,:))         CALL DEBUG_STATS_RL(Nr,GvNm1,'GvNm1 (DYNAMICS)',myThid)
814  cblk  write(0,*) 'dynamics: K23',minval(K23(1:sNx,1:sNy,:)),         CALL DEBUG_STATS_RL(Nr,GtNm1,'GtNm1 (DYNAMICS)',myThid)
815  cblk &                           maxval(K23(1:sNx,1:sNy,:))         CALL DEBUG_STATS_RL(Nr,GsNm1,'GsNm1 (DYNAMICS)',myThid)
816  cblk  write(0,*) 'dynamics: K33',minval(K33(1:sNx,1:sNy,:)),        ENDIF
817  cblk &                           maxval(K33(1:sNx,1:sNy,:))  #endif
 C     write(0,*) 'dynamics: gT ',minval(gT(1:sNx,1:sNy,:,:,:)),  
 C    &                           maxval(gT(1:sNx,1:sNy,:,:,:))  
 C     write(0,*) 'dynamics: T  ',minval(Theta(1:sNx,1:sNy,:,:,:)),  
 C    &                           maxval(Theta(1:sNx,1:sNy,:,:,:))  
 C     write(0,*) 'dynamics: gS ',minval(gS(1:sNx,1:sNy,:,:,:)),  
 C    &                           maxval(gS(1:sNx,1:sNy,:,:,:))  
 C     write(0,*) 'dynamics: S  ',minval(salt(1:sNx,1:sNy,:,:,:)),  
 C    &                           maxval(salt(1:sNx,1:sNy,:,:,:))  
 C     write(0,*) 'dynamics: pH ',minval(pH/(Gravity*Rhonil),mask=ph.NE.0.),  
 C    &                           maxval(pH/(Gravity*Rhonil))  
818    
819        RETURN        RETURN
820        END        END
Legend:
Removed from v.1.20  
changed lines
  Added in v.1.74
  ViewVC Help
Powered by ViewVC 1.1.22