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