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