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