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