/[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.49 by heimbach, Fri Jun 9 02:45:04 2000 UTC revision 1.65 by jmc, Thu Mar 8 20:25:01 2001 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2    C $Name$
3    
4  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
5    
# 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     \==========================================================/
 c  
 c     changed: Patrick Heimbach heimbach@mit.edu 6-Jun-2000  
 c              - computation of ikey wrong for nTx,nTy > 1  
 c                and/or nsx,nsy > 1: act1 and act2 were  
 c                mixed up.  
   
24        IMPLICIT NONE        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"  #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  #ifdef ALLOW_KPP
39  #include "KPPMIX.h"  # include "KPP.h"
40  #endif  #endif
41    
42  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_TIMEAVE
43  #include "tamc.h"  #include "TIMEAVE_STATV.h"
 #include "tamc_keys.h"  
44  #endif  #endif
45    
46  C     == Routine arguments ==  C     == Routine arguments ==
# Line 57  C     == Local variables Line 55  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, rTrans - Per block temporaries holding flow  C     uTrans, vTrans, rTrans - Per block temporaries holding flow
57  C                              transport  C                              transport
58  C     rVel                     o uTrans: Zonal transport  C                              o uTrans: Zonal transport
59  C                              o vTrans: Meridional transport  C                              o vTrans: Meridional transport
60  C                              o rTrans: Vertical transport  C                              o rTrans: Vertical transport
 C                              o rVel:   Vertical velocity at upper and  
 C                                        lower 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     rhoK, rhoKM1   - Density at current level, level above and level  C     rhoK, rhoKM1   - Density at current level, and level above
 C                      below.  
 C     rhoKP1                                                                    
 C     buoyK, buoyKM1 - Buoyancy at current level and level above.  
68  C     phiHyd         - Hydrostatic part of the potential phiHydi.  C     phiHyd         - Hydrostatic part of the potential phiHydi.
69  C                      In z coords phiHydiHyd is the hydrostatic  C                      In z coords phiHydiHyd is the hydrostatic
70  C                      pressure anomaly  C                      Potential (=pressure/rho0) anomaly
71  C                      In p coords phiHydiHyd is the geopotential  C                      In p coords phiHydiHyd is the geopotential
72  C                      surface height  C                      surface height anomaly.
73  C                      anomaly.  C     phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean)
74  C     etaSurfX,      - Holds surface elevation gradient in X and Y.  C     phiSurfY             or geopotentiel (atmos) in X and Y direction
 C     etaSurfY  
 C     K13, K23, K33  - Non-zero elements of small-angle approximation  
 C                      diffusion tensor.  
 C     KapGM          - Spatially varying Visbeck et. al mixing coeff.  
75  C     KappaRT,       - Total diffusion in vertical for T and S.  C     KappaRT,       - Total diffusion in vertical for T and S.
76  C     KappaRS          (background + spatially varying, isopycnal term).  C     KappaRS          (background + spatially varying, isopycnal term).
77  C     iMin, iMax     - Ranges and sub-block indices on which calculations  C     iMin, iMax     - Ranges and sub-block indices on which calculations
78  C     jMin, jMax       are applied.  C     jMin, jMax       are applied.
79  C     bi, bj  C     bi, bj
80  C     k, kUp,        - Index for layer above and below. kUp and kDown  C     k, kup,        - Index for layer above and below. kup and kDown
81  C     kDown, kM1       are switched with layer to be the appropriate  C     kDown, km1       are switched with layer to be the appropriate
82  C                      index 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 rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
       _RL rVel    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)  
88        _RS maskC   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS maskC   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
89        _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
       _RL aTerm   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL xTerm   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL cTerm   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL mTerm   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL pTerm   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL fZon    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL fMer    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
90        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
91        _RL fVerS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
92        _RL fVerU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
93        _RL fVerV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
94        _RL phiHyd  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL phiHyd  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
95        _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
       _RL rhokp1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
96        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
97        _RL buoyKM1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
98        _RL buoyK   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
       _RL rhotmp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL etaSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL etaSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL K13     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
       _RL K23     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
       _RL K33     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
       _RL KapGM   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
99        _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
100        _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
101        _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
102        _RL KappaRV (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRV (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
103          _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
104          _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
105          _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
106    
107  #ifdef INCLUDE_CONVECT_CALL  C This is currently used by IVDC and Diagnostics
108        _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
 #endif  
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  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
123        INTEGER    isbyte        INTEGER    isbyte
124        PARAMETER( isbyte = 4 )        PARAMETER( isbyte = 4 )
125    
126        INTEGER act1, act2, act3, act4        INTEGER act1, act2, act3, act4
127        INTEGER max1, max2, max3        INTEGER max1, max2, max3
128        INTEGER ikact, iikey,kkey        INTEGER iikey, kkey
129        INTEGER maximpl        INTEGER maximpl
130  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
131    
132  C---    The algorithm...  C---    The algorithm...
133  C  C
# Line 171  C Line 142  C
142  C       "Calculation of Gs"  C       "Calculation of Gs"
143  C       ===================  C       ===================
144  C       This is where all the accelerations and tendencies (ie.  C       This is where all the accelerations and tendencies (ie.
145  C       phiHydysics, parameterizations etc...) are calculated  C       physics, parameterizations etc...) are calculated
 C         rVel = sum_r ( div. u[n] )  
146  C         rho = rho ( theta[n], salt[n] )  C         rho = rho ( theta[n], salt[n] )
147  C         b   = b(rho, theta)  C         b   = b(rho, theta)
148  C         K31 = K31 ( rho )  C         K31 = K31 ( rho )
149  C         Gu[n] = Gu( u[n], v[n], rVel, b, ... )  C         Gu[n] = Gu( u[n], v[n], wVel, b, ... )
150  C         Gv[n] = Gv( u[n], v[n], rVel, b, ... )  C         Gv[n] = Gv( u[n], v[n], wVel, b, ... )
151  C         Gt[n] = Gt( theta[n], u[n], v[n], rVel, K31, ... )  C         Gt[n] = Gt( theta[n], u[n], v[n], wVel, K31, ... )
152  C         Gs[n] = Gs( salt[n], u[n], v[n], rVel, K31, ... )  C         Gs[n] = Gs( salt[n], u[n], v[n], wVel, K31, ... )
153  C  C
154  C       "Time-stepping" or "Prediction"  C       "Time-stepping" or "Prediction"
155  C       ================================  C       ================================
# Line 206  C--- Line 176  C---
176  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
177  C--   dummy statement to end declaration part  C--   dummy statement to end declaration part
178        ikey = 1        ikey = 1
179  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
180    
181  C--   Set up work arrays with valid (i.e. not NaN) values  C--   Set up work arrays with valid (i.e. not NaN) values
182  C     These inital values do not alter the numerical results. They  C     These inital values do not alter the numerical results. They
# Line 219  C     uninitialised but inert locations. Line 189  C     uninitialised but inert locations.
189          yA(i,j)      = 0. _d 0          yA(i,j)      = 0. _d 0
190          uTrans(i,j)  = 0. _d 0          uTrans(i,j)  = 0. _d 0
191          vTrans(i,j)  = 0. _d 0          vTrans(i,j)  = 0. _d 0
192          aTerm(i,j)   = 0. _d 0          DO k=1,Nr
193          xTerm(i,j)   = 0. _d 0           phiHyd(i,j,k)  = 0. _d 0
         cTerm(i,j)   = 0. _d 0  
         mTerm(i,j)   = 0. _d 0  
         pTerm(i,j)   = 0. _d 0  
         fZon(i,j)    = 0. _d 0  
         fMer(i,j)    = 0. _d 0  
         DO K=1,Nr  
          phiHyd (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  
194           KappaRU(i,j,k) = 0. _d 0           KappaRU(i,j,k) = 0. _d 0
195           KappaRV(i,j,k) = 0. _d 0           KappaRV(i,j,k) = 0. _d 0
196             sigmaX(i,j,k) = 0. _d 0
197             sigmaY(i,j,k) = 0. _d 0
198             sigmaR(i,j,k) = 0. _d 0
199          ENDDO          ENDDO
200          rhoKM1 (i,j) = 0. _d 0          rhoKM1 (i,j) = 0. _d 0
201          rhok   (i,j) = 0. _d 0          rhok   (i,j) = 0. _d 0
         rhoKP1 (i,j) = 0. _d 0  
         rhoTMP (i,j) = 0. _d 0  
         buoyKM1(i,j) = 0. _d 0  
         buoyK  (i,j) = 0. _d 0  
202          maskC  (i,j) = 0. _d 0          maskC  (i,j) = 0. _d 0
203            phiSurfX(i,j) = 0. _d 0
204            phiSurfY(i,j) = 0. _d 0
205         ENDDO         ENDDO
206        ENDDO        ENDDO
207    
208    
209  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
210  C--   HPF directive to help TAMC  C--   HPF directive to help TAMC
211  !HPF$ INDEPENDENT  CHPF$ INDEPENDENT
212  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
213    
214        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
215    
216  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
217  C--    HPF directive to help TAMC  C--    HPF directive to help TAMC
218  !HPF$  INDEPENDENT, NEW (rTrans,rVel,fVerT,fVerS,fVerU,fVerV  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS,fVerU,fVerV
219  !HPF$&                  ,phiHyd,K13,K23,K33,KapGM  CHPF$&                  ,phiHyd,utrans,vtrans,maskc,xA,yA
220  !HPF$&                  ,utrans,vtrans,maskc,xA,yA  CHPF$&                  ,KappaRT,KappaRS,KappaRU,KappaRV
221  !HPF$&                  ,KappaRT,KappaRS,KappaRU,KappaRV  CHPF$&                  )
222  !HPF$&                  )  #endif /* ALLOW_AUTODIFF_TAMC */
 #endif  
223    
224         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
225    
# Line 278  C--    HPF directive to help TAMC Line 238  C--    HPF directive to help TAMC
238            ikey = (act1 + 1) + act2*max1            ikey = (act1 + 1) + act2*max1
239       &                      + act3*max1*max2       &                      + act3*max1*max2
240       &                      + act4*max1*max2*max3       &                      + act4*max1*max2*max3
241  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
242    
243  C--     Set up work arrays that need valid initial values  C--     Set up work arrays that need valid initial values
244          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
245           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
246            rTrans(i,j)   = 0. _d 0            rTrans(i,j)   = 0. _d 0
           rVel  (i,j,1) = 0. _d 0  
           rVel  (i,j,2) = 0. _d 0  
247            fVerT (i,j,1) = 0. _d 0            fVerT (i,j,1) = 0. _d 0
248            fVerT (i,j,2) = 0. _d 0            fVerT (i,j,2) = 0. _d 0
249            fVerS (i,j,1) = 0. _d 0            fVerS (i,j,1) = 0. _d 0
# Line 294  C--     Set up work arrays that need val Line 252  C--     Set up work arrays that need val
252            fVerU (i,j,2) = 0. _d 0            fVerU (i,j,2) = 0. _d 0
253            fVerV (i,j,1) = 0. _d 0            fVerV (i,j,1) = 0. _d 0
254            fVerV (i,j,2) = 0. _d 0            fVerV (i,j,2) = 0. _d 0
           phiHyd(i,j,1) = 0. _d 0  
           K13   (i,j,1) = 0. _d 0  
           K23   (i,j,1) = 0. _d 0  
           K33   (i,j,1) = 0. _d 0  
           KapGM (i,j)   = GMkbackground  
255           ENDDO           ENDDO
256          ENDDO          ENDDO
257    
258          DO k=1,Nr          DO k=1,Nr
259           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
260            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
261  #ifdef INCLUDE_CONVECT_CALL  C This is currently also used by IVDC and Diagnostics
262             ConvectCount(i,j,k) = 0.             ConvectCount(i,j,k) = 0.
 #endif  
263             KappaRT(i,j,k) = 0. _d 0             KappaRT(i,j,k) = 0. _d 0
264             KappaRS(i,j,k) = 0. _d 0             KappaRS(i,j,k) = 0. _d 0
265            ENDDO            ENDDO
# Line 320  C--     Set up work arrays that need val Line 272  C--     Set up work arrays that need val
272          jMax = sNy+OLy          jMax = sNy+OLy
273    
274    
275          K = 1  C--     Start of diagnostic loop
276          BOTTOM_LAYER = K .EQ. Nr          DO k=Nr,1,-1
   
 #ifdef DO_PIPELINED_CORRECTION_STEP  
 C--     Calculate gradient of surface pressure  
         CALL CALC_GRAD_ETA_SURF(  
      I       bi,bj,iMin,iMax,jMin,jMax,  
      O       etaSurfX,etaSurfY,  
      I       myThid)  
 C--     Update fields in top level according to tendency terms  
         CALL CORRECTION_STEP(  
      I       bi,bj,iMin,iMax,jMin,jMax,K,  
      I       etaSurfX,etaSurfY,myTime,myThid)  
   
 #ifdef ALLOW_OBCS  
         IF (openBoundaries) THEN  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE uvel (:,:,k,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  
 CADJ STORE vvel (:,:,k,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  
 CADJ STORE theta(:,:,k,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  
 CADJ STORE salt(:,:,k,bi,bj)   = comlev1_2d, key = ikey, byte = isbyte  
 #endif  
            CALL APPLY_OBCS1( bi, bj, K, myThid )  
         END IF  
 #endif  
277    
         IF ( .NOT. BOTTOM_LAYER ) THEN  
 C--      Update fields in layer below according to tendency terms  
          CALL CORRECTION_STEP(  
      I        bi,bj,iMin,iMax,jMin,jMax,K+1,  
      I        etaSurfX,etaSurfY,myTime,myThid)  
 #ifdef ALLOW_OBCS  
          IF (openBoundaries) THEN  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE uvel (:,:,k,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  
 CADJ STORE vvel (:,:,k,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  
 CADJ STORE theta(:,:,k,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  
 CADJ STORE salt(:,:,k,bi,bj)   = comlev1_2d, key = ikey, byte = isbyte  
 #endif  
             CALL APPLY_OBCS1( bi, bj, K+1, myThid )  
          END IF  
 #endif  
         ENDIF  
 #endif  
 C--     Density of 1st level (below W(1)) reference to level 1  
 #ifdef  INCLUDE_FIND_RHO_CALL  
278  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
279  CADJ STORE theta(:,:,k,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  C? Patrick, is this formula correct now that we change the loop range?
280  CADJ STORE salt (:,:,k,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  C? Do we still need this?
281  #endif           kkey = (ikey-1)*(Nr-2+1) + (k-2) + 1
282          CALL FIND_RHO(  #endif /* ALLOW_AUTODIFF_TAMC */
283       I     bi, bj, iMin, iMax, jMin, jMax, K, K, eosType,  
284       O     rhoKm1,  C--       Integrate continuity vertically for vertical velocity
285       I     myThid )            CALL INTEGRATE_FOR_W(
286  #endif       I                         bi, bj, k, uVel, vVel,
287         O                         wVel,
288          IF (       (.NOT. BOTTOM_LAYER)       I                         myThid )
289  #ifdef ALLOW_KPP  
290       &       .AND. (.NOT.usingKPPmixing) ! CONVECT not needed with KPP mixing  #ifdef    ALLOW_OBCS
291  #endif  #ifdef    ALLOW_NONHYDROSTATIC
292       &     ) THEN  C--       Apply OBC to W if in N-H mode
293  C--      Check static stability with layer below            IF (useOBCS.AND.nonHydrostatic) THEN
294  C--      and mix as needed.              CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )
295  #ifdef  INCLUDE_FIND_RHO_CALL            ENDIF
296  #ifdef ALLOW_AUTODIFF_TAMC  #endif    /* ALLOW_NONHYDROSTATIC */
297  CADJ STORE theta(:,:,k,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  #endif    /* ALLOW_OBCS */
298  CADJ STORE salt (:,:,k,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  
299  #endif  C--       Calculate gradients of potential density for isoneutral
300           CALL FIND_RHO(  C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)
301       I      bi, bj, iMin, iMax, jMin, jMax, K+1, K, eosType,  c         IF ( k.GT.1 .AND. (useGMRedi.OR.ivdc_kappa.NE.0.) ) THEN
302       O      rhoKp1,            IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.) ) THEN
303       I      myThid )              CALL FIND_RHO(
304  #endif       I        bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,
305         I        theta, salt,
306  #ifdef  INCLUDE_CONVECT_CALL       O        rhoK,
307         I        myThid )
308  #ifdef ALLOW_AUTODIFF_TAMC              IF (k.GT.1) CALL FIND_RHO(
309  CADJ STORE rhoKm1(:,:)  = comlev1_2d, key = ikey, byte = isbyte       I        bi, bj, iMin, iMax, jMin, jMax, k-1, k, eosType,
310  CADJ STORE rhoKp1(:,:)  = comlev1_2d, key = ikey, byte = isbyte       I        theta, salt,
311  #endif       O        rhoKm1,
312           CALL CONVECT(       I        myThid )
313       I       bi,bj,iMin,iMax,jMin,jMax,K+1,rhoKm1,rhoKp1,              CALL GRAD_SIGMA(
314       U       ConvectCount,       I             bi, bj, iMin, iMax, jMin, jMax, k,
315       I       myTime,myIter,myThid)       I             rhoK, rhoKm1, rhoK,
316  #ifdef ALLOW_AUTODIFF_TAMC       O             sigmaX, sigmaY, sigmaR,
317  CADJ STORE theta(:,:,k+1,bi,bj),theta(:,:,k,bi,bj)       I             myThid )
318  CADJ &     = comlev1_2d, key = ikey, byte = isbyte            ENDIF
319  CADJ STORE salt (:,:,k+1,bi,bj),salt (:,:,k,bi,bj)  
320  CADJ &     = comlev1_2d, key = ikey, byte = isbyte  C--       Implicit Vertical Diffusion for Convection
321  #endif  c ==> should use sigmaR !!!
322              IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN
323                CALL CALC_IVDC(
324         I        bi, bj, iMin, iMax, jMin, jMax, k,
325         I        rhoKm1, rhoK,
326         U        ConvectCount, KappaRT, KappaRS,
327         I        myTime, myIter, myThid)
328              ENDIF
329    
330  #endif  C--     end of diagnostic k loop (Nr:1)
331            ENDDO
332    
333  C--      Implicit Vertical Diffusion for Convection  #ifdef  ALLOW_OBCS
334           IF (ivdc_kappa.NE.0.) CALL CALC_IVDC(  C--     Calculate future values on open boundaries
335       I       bi,bj,iMin,iMax,jMin,jMax,K+1,rhoKm1,rhoKp1,          IF (useOBCS) THEN
336       U       ConvectCount, KappaRT, KappaRS,            CALL OBCS_CALC( bi, bj, myTime+deltaT,
337       I       myTime,myIter,myThid)       I            uVel, vVel, wVel, theta, salt,
338  CRG: do we need do store STORE KappaRT, KappaRS ?       I            myThid )
   
 C--      Recompute density after mixing  
 #ifdef  INCLUDE_FIND_RHO_CALL  
          CALL FIND_RHO(  
      I      bi, bj, iMin, iMax, jMin, jMax, K, K, eosType,  
      O      rhoKm1,  
      I      myThid )  
 #endif  
339          ENDIF          ENDIF
340  C--     Calculate buoyancy  #endif  /* ALLOW_OBCS */
         CALL CALC_BUOYANCY(  
      I      bi,bj,iMin,iMax,jMin,jMax,K,rhoKm1,  
      O      buoyKm1,  
      I      myThid )  
 C--     Integrate hydrostatic balance for phiHyd with BC of  
 C--     phiHyd(z=0)=0  
         CALL CALC_PHI_HYD(  
      I      bi,bj,iMin,iMax,jMin,jMax,K,buoyKm1,buoyKm1,  
      U      phiHyd,  
      I      myThid )  
   
 C----------------------------------------------  
 C--     start of downward loop  
 C----------------------------------------------  
         DO K=2,Nr  
341    
342    C--     Determines forcing terms based on external fields
343    C       relaxation terms, etc.
344            CALL EXTERNAL_FORCING_SURF(
345         I             bi, bj, iMin, iMax, jMin, jMax,
346         I             myThid )
347    
348    #ifdef  ALLOW_GMREDI
349    C--     Calculate iso-neutral slopes for the GM/Redi parameterisation
350            IF (useGMRedi) THEN
351              DO k=1,Nr
352                CALL GMREDI_CALC_TENSOR(
353         I             bi, bj, iMin, iMax, jMin, jMax, k,
354         I             sigmaX, sigmaY, sigmaR,
355         I             myThid )
356              ENDDO
357  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
358           kkey = ikact*(Nr-2+1) + (k-2) + 1          ELSE
359  #endif            DO k=1, Nr
360                CALL GMREDI_CALC_TENSOR_DUMMY(
361           BOTTOM_LAYER = K .EQ. Nr       I             bi, bj, iMin, iMax, jMin, jMax, k,
362         I             sigmaX, sigmaY, sigmaR,
363  #ifdef DO_PIPELINED_CORRECTION_STEP       I             myThid )
364           IF ( .NOT. BOTTOM_LAYER ) THEN            ENDDO
365  C--       Update fields in layer below according to tendency terms  #endif /* ALLOW_AUTODIFF_TAMC */
366            CALL CORRECTION_STEP(          ENDIF
367       I         bi,bj,iMin,iMax,jMin,jMax,K+1,  #endif  /* ALLOW_GMREDI */
      I         etaSurfX,etaSurfY,myTime,myThid)  
 #ifdef ALLOW_OBCS  
           IF (openBoundaries) THEN  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE uvel (:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  
 CADJ STORE vvel (:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  
 CADJ STORE theta(:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  
 CADJ STORE salt(:,:,k,bi,bj)   = comlev1_2d, key = ikey, byte = isbyte  
 #endif  
              CALL APPLY_OBCS1( bi, bj, K+1, myThid )  
           END IF  
 #endif  
          ENDIF  
 #endif  
368    
369  C--      Density of K level (below W(K)) reference to K level  #ifdef  ALLOW_KPP
370  #ifdef  INCLUDE_FIND_RHO_CALL  C--     Compute KPP mixing coefficients
371  #ifdef ALLOW_AUTODIFF_TAMC          IF (useKPP) THEN
372  CADJ STORE theta(:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte            CALL KPP_CALC(
373  CADJ STORE salt (:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte       I                  bi, bj, myTime, myThid )
374  #endif          ENDIF
375           CALL FIND_RHO(  #endif  /* ALLOW_KPP */
      I      bi, bj, iMin, iMax, jMin, jMax,  K, K, eosType,  
      O      rhoK,  
      I      myThid )  
 #endif  
          IF (       (.NOT. BOTTOM_LAYER)  
 #ifdef ALLOW_KPP  
      &       .AND. (.NOT.usingKPPmixing) ! CONVECT not needed with KPP mixing  
 #endif  
      &      ) THEN  
 C--       Check static stability with layer below and mix as needed.  
 C--       Density of K+1 level (below W(K+1)) reference to K level.  
 #ifdef  INCLUDE_FIND_RHO_CALL  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE theta(:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  
 CADJ STORE salt (:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  
 #endif  
           CALL FIND_RHO(  
      I       bi, bj, iMin, iMax, jMin, jMax,  K+1, K, eosType,  
      O       rhoKp1,  
      I       myThid )  
 #endif  
376    
377  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
378  CADJ STORE rhok  (:,:)   = comlev1_3d, key = kkey, byte = isbyte  CADJ STORE KappaRT(:,:,:)     = comlev1_bibj, key = ikey, byte = isbyte
379  CADJ STORE rhoKm1(:,:)   = comlev1_3d, key = kkey, byte = isbyte  CADJ STORE KappaRS(:,:,:)     = comlev1_bibj, key = ikey, byte = isbyte
380  CADJ STORE rhoKp1(:,:)   = comlev1_3d, key = kkey, byte = isbyte  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
381  #endif  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
382    CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
383  #ifdef  INCLUDE_CONVECT_CALL  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
384            CALL CONVECT(  #endif /* ALLOW_AUTODIFF_TAMC */
385       I        bi,bj,iMin,iMax,jMin,jMax,K+1,rhoK,rhoKp1,  
386       U        ConvectCount,  #ifdef ALLOW_AIM
387       I        myTime,myIter,myThid)  C       AIM - atmospheric intermediate model, physics package code.
388  #ifdef ALLOW_AUTODIFF_TAMC  C note(jmc) : phiHyd=0 at this point but is not really used in Molteni Physics
389  CADJ STORE theta(:,:,k+1,bi,bj),theta(:,:,k,bi,bj)          IF ( useAIM ) THEN
390  CADJ &     = comlev1_3d, key = kkey, byte = isbyte           CALL TIMER_START('AIM_DO_ATMOS_PHYS      [DYNAMICS]', myThid)
391  CADJ STORE salt (:,:,k+1,bi,bj),salt (:,:,k,bi,bj)           CALL AIM_DO_ATMOS_PHYSICS( phiHyd, myTime, myThid )
392  CADJ &     = comlev1_3d, key = kkey, byte = isbyte           CALL TIMER_STOP ('AIM_DO_ATMOS_PHYS      [DYNAMICS]', myThid)
393  #endif          ENDIF
394  #endif  #endif /* ALLOW_AIM */
   
 C--      Implicit Vertical Diffusion for Convection  
          IF (ivdc_kappa.NE.0.) THEN  
             CALL CALC_IVDC(  
      I       bi,bj,iMin,iMax,jMin,jMax,K+1,rhoKm1,rhoKp1,  
      U       ConvectCount, KappaRT, KappaRS,  
      I       myTime,myIter,myThid)  
 CRG: do we need do store STORE KappaRT, KappaRS ?  
          END IF  
   
 C--       Recompute density after mixing  
 #ifdef  INCLUDE_FIND_RHO_CALL  
           CALL FIND_RHO(  
      I       bi, bj, iMin, iMax, jMin, jMax, K, K, eosType,  
      O       rhoK,  
      I       myThid )  
 #endif  
          ENDIF  
 C--      Calculate buoyancy  
          CALL CALC_BUOYANCY(  
      I       bi,bj,iMin,iMax,jMin,jMax,K,rhoK,  
      O       buoyK,  
      I       myThid )  
 C--      Integrate hydrostatic balance for phiHyd with BC of  
 C--      phiHyd(z=0)=0  
          CALL CALC_PHI_HYD(  
      I        bi,bj,iMin,iMax,jMin,jMax,K,buoyKm1,buoyK,  
      U        phiHyd,  
      I        myThid )  
 C--      Calculate iso-neutral slopes for the GM/Redi parameterisation  
 #ifdef  INCLUDE_FIND_RHO_CALL  
          CALL FIND_RHO(  
      I        bi, bj, iMin, iMax, jMin, jMax, K-1, K, eosType,  
      O        rhoTmp,  
      I        myThid )  
 #endif  
395    
 #ifdef  INCLUDE_CALC_ISOSLOPES_CALL  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE rhoTmp(:,:)  = comlev1_3d, key = kkey, byte = isbyte  
 CADJ STORE rhok  (:,:)  = comlev1_3d, key = kkey, byte = isbyte  
 CADJ STORE rhoKm1(:,:)  = comlev1_3d, key = kkey, byte = isbyte  
 CADJ STORE kapgm (:,:)  = comlev1_3d, key = kkey, byte = isbyte  
 #endif  
          CALL CALC_ISOSLOPES(  
      I        bi, bj, iMin, iMax, jMin, jMax, K,  
      I        rhoKm1, rhoK, rhotmp,  
      O        K13, K23, K33, KapGM,  
      I        myThid )  
 #endif  
396    
397           DO J=jMin,jMax  C--     Start of thermodynamics loop
398            DO I=iMin,iMax          DO k=Nr,1,-1
 #ifdef  INCLUDE_FIND_RHO_CALL  
            rhoKm1 (I,J) = rhoK(I,J)  
 #endif  
            buoyKm1(I,J) = buoyK(I,J)  
           ENDDO  
          ENDDO  
         ENDDO  
 C--     end of k loop  
399    
400  #ifdef ALLOW_AUTODIFF_TAMC  C--       km1    Points to level above k (=k-1)
401  CADJ STORE theta(:,:,:,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  C--       kup    Cycles through 1,2 to point to layer above
402  CADJ STORE salt (:,:,:,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  C--       kDown  Cycles through 2,1 to point to current layer
 CADJ STORE uvel (:,:,:,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  
 CADJ STORE vvel (:,:,:,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  
 #endif  
403    
404  #ifdef ALLOW_KPP            km1  = MAX(1,k-1)
405  C----------------------------------------------            kup  = 1+MOD(k+1,2)
406  C--     Compute KPP mixing coefficients            kDown= 1+MOD(k,2)
 C----------------------------------------------  
         IF (usingKPPmixing) THEN  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE fu  (:,:  ,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  
 CADJ STORE fv  (:,:  ,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  
 #endif  
          CALL TIMER_START('KVMIX (FIND KPP COEFFICIENTS) [DYNAMICS]'  
      I          , myThid)  
          CALL KVMIX(  
      I               bi, bj, myTime, myThid )  
          CALL TIMER_STOP ('KVMIX (FIND KPP COEFFICIENTS) [DYNAMICS]'  
      I        , myThid)  
         ENDIF  
 #endif  
407    
408  C----------------------------------------------            iMin = 1-OLx+2
409  C--     start of upward loop            iMax = sNx+OLx-1
410  C----------------------------------------------            jMin = 1-OLy+2
411          DO K = Nr, 1, -1            jMax = sNy+OLy-1
   
          kM1  =max(1,k-1)   ! Points to level above k (=k-1)  
          kUp  =1+MOD(k+1,2) ! Cycles through 1,2 to point to layer above  
          kDown=1+MOD(k,2)   ! Cycles through 2,1 to point to current layer  
   
          iMin = 1-OLx+2  
          iMax = sNx+OLx-1  
          jMin = 1-OLy+2  
          jMax = sNy+OLy-1  
412    
413  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
414           kkey = ikact*(Nr-1+1) + (k-1) + 1  CPatrick Is this formula correct?
415  #endif           kkey = (ikey-1)*(Nr-1+1) + (k-1) + 1
416    CADJ STORE rTrans(:,:)       = comlev1_bibj_k, key = kkey, byte = isbyte
417  #ifdef ALLOW_AUTODIFF_TAMC  CADJ STORE KappaRT(:,:,:)    = comlev1_bibj_k, key = kkey, byte = isbyte
418  CADJ STORE rvel  (:,:,kDown)  = comlev1_3d, key = kkey, byte = isbyte  CADJ STORE KappaRS(:,:,:)    = comlev1_bibj_k, key = kkey, byte = isbyte
419  CADJ STORE rTrans(:,:)        = comlev1_3d, key = kkey, byte = isbyte  #endif /* ALLOW_AUTODIFF_TAMC */
 CADJ STORE KappaRT(:,:,:)     = comlev1_3d, key = kkey, byte = isbyte  
 CADJ STORE KappaRS(:,:,:)     = comlev1_3d, key = kkey, byte = isbyte  
 #endif  
420    
421  C--      Get temporary terms used by tendency routines  C--      Get temporary terms used by tendency routines
422           CALL CALC_COMMON_FACTORS (           CALL CALC_COMMON_FACTORS (
423       I        bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,       I        bi,bj,iMin,iMax,jMin,jMax,k,km1,kup,kDown,
424       O        xA,yA,uTrans,vTrans,rTrans,rVel,maskC,maskUp,       O        xA,yA,uTrans,vTrans,rTrans,maskC,maskUp,
425       I        myThid)       I        myThid)
426    
 #ifdef ALLOW_OBCS  
         IF (openBoundaries) THEN  
          CALL APPLY_OBCS3( bi, bj, K, Kup, rTrans, rVel, myThid )  
         ENDIF  
 #endif  
   
427  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL
428  C--      Calculate the total vertical diffusivity  C--      Calculate the total vertical diffusivity
429           CALL CALC_DIFFUSIVITY(           CALL CALC_DIFFUSIVITY(
430       I        bi,bj,iMin,iMax,jMin,jMax,K,       I        bi,bj,iMin,iMax,jMin,jMax,k,
431       I        maskC,maskUp,KapGM,K33,       I        maskC,maskup,
432       O        KappaRT,KappaRS,KappaRU,KappaRV,       O        KappaRT,KappaRS,KappaRU,KappaRV,
433       I        myThid)       I        myThid)
434  #endif  #endif
435  C--      Calculate accelerations in the momentum equations  
436           IF ( momStepping ) THEN  C--      Calculate active tracer tendencies (gT,gS,...)
437            CALL CALC_MOM_RHS(  C        and step forward storing result in gTnm1, gSnm1, etc.
      I         bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,  
      I         xA,yA,uTrans,vTrans,rTrans,rVel,maskC,  
      I         phiHyd,KappaRU,KappaRV,  
      U         aTerm,xTerm,cTerm,mTerm,pTerm,  
      U         fZon, fMer, fVerU, fVerV,  
      I         myTime, myThid)  
 #ifdef ALLOW_AUTODIFF_TAMC  
 #ifdef INCLUDE_CD_CODE  
          ELSE  
             DO j=1-OLy,sNy+OLy  
                DO i=1-OLx,sNx+OLx  
                   guCD(i,j,k,bi,bj) = 0.0  
                   gvCD(i,j,k,bi,bj) = 0.0  
                END DO  
             END DO  
 #endif  
 #endif  
          ENDIF  
 C--      Calculate active tracer tendencies  
438           IF ( tempStepping ) THEN           IF ( tempStepping ) THEN
439            CALL CALC_GT(             CALL CALC_GT(
440       I         bi,bj,iMin,iMax,jMin,jMax, k,kM1,kUp,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
441       I         xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,       I         xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,
442       I         K13,K23,KappaRT,KapGM,       I         KappaRT,
443       U         aTerm,xTerm,fZon,fMer,fVerT,       U         fVerT,
444       I         myTime, myThid)       I         myTime, myThid)
445               CALL TIMESTEP_TRACER(
446         I         bi,bj,iMin,iMax,jMin,jMax,k,
447         I         theta, gT,
448         U         gTnm1,
449         I         myIter, myThid)
450           ENDIF           ENDIF
451           IF ( saltStepping ) THEN           IF ( saltStepping ) THEN
452            CALL CALC_GS(             CALL CALC_GS(
453       I         bi,bj,iMin,iMax,jMin,jMax, k,kM1,kUp,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
454       I         xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,       I         xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,
455       I         K13,K23,KappaRS,KapGM,       I         KappaRS,
456       U         aTerm,xTerm,fZon,fMer,fVerS,       U         fVerS,
457       I         myTime, myThid)       I         myTime, myThid)
458               CALL TIMESTEP_TRACER(
459         I         bi,bj,iMin,iMax,jMin,jMax,k,
460         I         salt, gS,
461         U         gSnm1,
462         I         myIter, myThid)
463           ENDIF           ENDIF
464  #ifdef ALLOW_OBCS  
465  C--      Calculate future values on open boundaries  #ifdef   ALLOW_OBCS
          IF (openBoundaries) THEN  
 Caja      CALL CYCLE_OBCS( K, bi, bj, myThid )  
           CALL SET_OBCS( K, bi, bj, myTime+deltaTclock, myThid )  
          ENDIF  
 #endif  
 C--      Prediction step (step forward all model variables)  
          CALL TIMESTEP(  
      I       bi,bj,iMin,iMax,jMin,jMax,K,  
      I       myIter, myThid)  
 #ifdef ALLOW_OBCS  
466  C--      Apply open boundary conditions  C--      Apply open boundary conditions
467           IF (openBoundaries) THEN           IF (useOBCS) THEN
468  #ifdef ALLOW_AUTODIFF_TAMC             CALL OBCS_APPLY_TS( bi, bj, k, gTnm1, gSnm1, myThid )
 CADJ STORE gunm1(:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  
 CADJ STORE gvnm1(:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  
 CADJ STORE gwnm1(:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  
 #endif  
             CALL APPLY_OBCS2( bi, bj, K, myThid )  
469           END IF           END IF
470  #endif  #endif   /* ALLOW_OBCS */
471    
472  C--      Freeze water  C--      Freeze water
473           IF (allowFreezing) THEN           IF (allowFreezing) THEN
474  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
475  CADJ STORE gTNm1(:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  CADJ STORE gTNm1(:,:,k,bi,bj) = comlev1_bibj_k
476  #endif  CADJ &   , key = kkey, byte = isbyte
477              CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, K, myThid )  #endif /* ALLOW_AUTODIFF_TAMC */
478                CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )
479           END IF           END IF
480    
481  #ifdef DIVG_IN_DYNAMICS  C--     end of thermodynamic k loop (Nr:1)
482  C--      Diagnose barotropic divergence of predicted fields          ENDDO
          CALL CALC_DIV_GHAT(  
      I       bi,bj,iMin,iMax,jMin,jMax,K,  
      I       xA,yA,  
      I       myThid)  
 #endif /* DIVG_IN_DYNAMICS */  
   
 C--      Cumulative diagnostic calculations (ie. time-averaging)  
 #ifdef INCLUDE_DIAGNOSTICS_INTERFACE_CODE  
          IF (taveFreq.GT.0.) THEN  
           CALL DO_TIME_AVERAGES(  
      I                           myTime, myIter, bi, bj, K, kUp, kDown,  
      I                           K13, K23, rVel, KapGM, ConvectCount,  
      I                           myThid )  
          ENDIF  
 #endif  
483    
484    
485          ENDDO ! K  #ifdef ALLOW_AUTODIFF_TAMC
486    CPatrick? What about this one?
487               maximpl = 6
488               iikey = (ikey-1)*maximpl
489    #endif /* ALLOW_AUTODIFF_TAMC */
490    
491  C--     Implicit diffusion  C--     Implicit diffusion
492          IF (implicitDiffusion) THEN          IF (implicitDiffusion) THEN
493    
 #ifdef ALLOW_AUTODIFF_TAMC  
            maximpl = 6  
            iikey = ikact*maximpl  
 #endif  
   
494           IF (tempStepping) THEN           IF (tempStepping) THEN
495  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
496              idkey = iikey + 1              idkey = iikey + 1
497  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
498              CALL IMPLDIFF(              CALL IMPLDIFF(
499       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
500       I         deltaTtracer, KappaRT,recip_HFacC,       I         deltaTtracer, KappaRT, recip_HFacC,
501       U         gTNm1,       U         gTNm1,
502       I         myThid )       I         myThid )
503           END IF           ENDIF
504    
505           IF (saltStepping) THEN           IF (saltStepping) THEN
506  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
507           idkey = iikey + 2           idkey = iikey + 2
508  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
509              CALL IMPLDIFF(              CALL IMPLDIFF(
510       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
511       I         deltaTtracer, KappaRS,recip_HFacC,       I         deltaTtracer, KappaRS, recip_HFacC,
512       U         gSNm1,       U         gSNm1,
513       I         myThid )       I         myThid )
514             ENDIF
515    
516    #ifdef   ALLOW_OBCS
517    C--      Apply open boundary conditions
518             IF (useOBCS) THEN
519               DO K=1,Nr
520                 CALL OBCS_APPLY_TS( bi, bj, k, gTnm1, gSnm1, myThid )
521               ENDDO
522           END IF           END IF
523    #endif   /* ALLOW_OBCS */
524    
525          ENDIF ! implicitDiffusion  C--     End If implicitDiffusion
526            ENDIF
527    
528  C--     Implicit viscosity  C--     Start computation of dynamics
529          IF (implicitViscosity) THEN          iMin = 1-OLx+2
530            iMax = sNx+OLx-1
531            jMin = 1-OLy+2
532            jMax = sNy+OLy-1
533    
534    C--     Explicit part of the Surface Potentiel Gradient (add in TIMESTEP)
535    C       (note: this loop will be replaced by CALL CALC_GRAD_ETA)
536            IF (implicSurfPress.NE.1.) THEN
537              CALL CALC_GRAD_PHI_SURF(
538         I         bi,bj,iMin,iMax,jMin,jMax,
539         I         etaN,
540         O         phiSurfX,phiSurfY,
541         I         myThid )                        
542            ENDIF
543    
544           IF (momStepping) THEN  C--     Start of dynamics loop
545  #ifdef ALLOW_AUTODIFF_TAMC          DO k=1,Nr
546           idkey = iikey + 3  
547  #endif  C--       km1    Points to level above k (=k-1)
548    C--       kup    Cycles through 1,2 to point to layer above
549    C--       kDown  Cycles through 2,1 to point to current layer
550    
551              km1  = MAX(1,k-1)
552              kup  = 1+MOD(k+1,2)
553              kDown= 1+MOD(k,2)
554    
555    C--      Integrate hydrostatic balance for phiHyd with BC of
556    C        phiHyd(z=0)=0
557    C        distinguishe between Stagger and Non Stagger time stepping
558             IF (staggerTimeStep) THEN
559               CALL CALC_PHI_HYD(
560         I        bi,bj,iMin,iMax,jMin,jMax,k,
561         I        gTnm1, gSnm1,
562         U        phiHyd,
563         I        myThid )
564             ELSE
565               CALL CALC_PHI_HYD(
566         I        bi,bj,iMin,iMax,jMin,jMax,k,
567         I        theta, salt,
568         U        phiHyd,
569         I        myThid )
570             ENDIF
571    
572    C--      Calculate accelerations in the momentum equations (gU, gV, ...)
573    C        and step forward storing the result in gUnm1, gVnm1, etc...
574             IF ( momStepping ) THEN
575               CALL CALC_MOM_RHS(
576         I         bi,bj,iMin,iMax,jMin,jMax,k,kup,kDown,
577         I         phiHyd,KappaRU,KappaRV,
578         U         fVerU, fVerV,
579         I         myTime, myThid)
580               CALL TIMESTEP(
581         I         bi,bj,iMin,iMax,jMin,jMax,k,
582         I         phiHyd, phiSurfX, phiSurfY,
583         I         myIter, myThid)
584    
585    #ifdef   ALLOW_OBCS
586    C--      Apply open boundary conditions
587             IF (useOBCS) THEN
588               CALL OBCS_APPLY_UV( bi, bj, k, gUnm1, gVnm1, myThid )
589             END IF
590    #endif   /* ALLOW_OBCS */
591    
592    #ifdef   ALLOW_AUTODIFF_TAMC
593    #ifdef   INCLUDE_CD_CODE
594             ELSE
595               DO j=1-OLy,sNy+OLy
596                 DO i=1-OLx,sNx+OLx
597                   guCD(i,j,k,bi,bj) = 0.0
598                   gvCD(i,j,k,bi,bj) = 0.0
599                 END DO
600               END DO
601    #endif   /* INCLUDE_CD_CODE */
602    #endif   /* ALLOW_AUTODIFF_TAMC */
603             ENDIF
604    
605    
606    C--     end of dynamics k loop (1:Nr)
607            ENDDO
608    
609    
610    
611    C--     Implicit viscosity
612            IF (implicitViscosity.AND.momStepping) THEN
613    #ifdef    ALLOW_AUTODIFF_TAMC
614              idkey = iikey + 3
615    #endif    /* ALLOW_AUTODIFF_TAMC */
616            CALL IMPLDIFF(            CALL IMPLDIFF(
617       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
618       I         deltaTmom, KappaRU,recip_HFacW,       I         deltaTmom, KappaRU,recip_HFacW,
619       U         gUNm1,       U         gUNm1,
620       I         myThid )       I         myThid )
621  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
622           idkey = iikey + 4            idkey = iikey + 4
623  #endif  #endif    /* ALLOW_AUTODIFF_TAMC */
624            CALL IMPLDIFF(            CALL IMPLDIFF(
625       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
626       I         deltaTmom, KappaRV,recip_HFacS,       I         deltaTmom, KappaRV,recip_HFacS,
627       U         gVNm1,       U         gVNm1,
628       I         myThid )       I         myThid )
629    
630  #ifdef INCLUDE_CD_CODE  #ifdef   ALLOW_OBCS
631    C--      Apply open boundary conditions
632             IF (useOBCS) THEN
633               DO K=1,Nr
634                 CALL OBCS_APPLY_UV( bi, bj, k, gUnm1, gVnm1, myThid )
635               ENDDO
636             END IF
637    #endif   /* ALLOW_OBCS */
638    
639  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef    INCLUDE_CD_CODE
640           idkey = iikey + 5  #ifdef    ALLOW_AUTODIFF_TAMC
641  #endif            idkey = iikey + 5
642    #endif    /* ALLOW_AUTODIFF_TAMC */
643            CALL IMPLDIFF(            CALL IMPLDIFF(
644       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
645       I         deltaTmom, KappaRU,recip_HFacW,       I         deltaTmom, KappaRU,recip_HFacW,
646       U         vVelD,       U         vVelD,
647       I         myThid )       I         myThid )
648  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef    ALLOW_AUTODIFF_TAMC
649          idkey = iikey + 6            idkey = iikey + 6
650  #endif  #endif    /* ALLOW_AUTODIFF_TAMC */
651            CALL IMPLDIFF(            CALL IMPLDIFF(
652       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
653       I         deltaTmom, KappaRV,recip_HFacS,       I         deltaTmom, KappaRV,recip_HFacS,
654       U         uVelD,       U         uVelD,
655       I         myThid )       I         myThid )
656    #endif    /* INCLUDE_CD_CODE */
657  #endif  C--     End If implicitViscosity.AND.momStepping
658            ENDIF
          ENDIF ! momStepping  
         ENDIF ! implicitViscosity  
659    
660    Cjmc : add for phiHyd output <- but not working if multi tile per CPU
661    c       IF ( DIFFERENT_MULTIPLE(dumpFreq,myTime+deltaTClock,myTime)
662    c    &  .AND. buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN
663    c         WRITE(suff,'(I10.10)') myIter+1
664    c         CALL WRITE_FLD_XYZ_RL('PH.',suff,phiHyd,myIter+1,myThid)
665    c       ENDIF
666    Cjmc(end)
667    
668    #ifdef ALLOW_TIMEAVE
669            IF (taveFreq.GT.0.) THEN
670              CALL TIMEAVE_CUMULATE(phiHydtave, phiHyd, Nr,
671         I                              deltaTclock, bi, bj, myThid)
672              IF (ivdc_kappa.NE.0.) THEN
673                CALL TIMEAVE_CUMULATE(ConvectCountTave, ConvectCount, Nr,
674         I                              deltaTclock, bi, bj, myThid)
675              ENDIF
676            ENDIF
677    #endif /* ALLOW_TIMEAVE */
678    
679         ENDDO         ENDDO
680        ENDDO        ENDDO
681    
 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: rVel(1) ',  
 C    &            minval(rVel(1:sNx,1:sNy,1),mask=rVel(1:sNx,1:sNy,1).NE.0.),  
 C    &            maxval(rVel(1:sNx,1:sNy,1),mask=rVel(1:sNx,1:sNy,1).NE.0.)  
 C     write(0,*) 'dynamics: rVel(2) ',  
 C    &            minval(rVel(1:sNx,1:sNy,2),mask=rVel(1:sNx,1:sNy,2).NE.0.),  
 C    &            maxval(rVel(1:sNx,1:sNy,2),mask=rVel(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: phiHyd ',minval(phiHyd/(Gravity*Rhonil),mask=phiHyd.NE.0.),  
 C    &                           maxval(phiHyd/(Gravity*Rhonil))  
 C     CALL PLOT_FIELD_XYZRL( gU, ' GU exiting dyanmics ' ,  
 C    &Nr, 1, myThid )  
 C     CALL PLOT_FIELD_XYZRL( gV, ' GV exiting dyanmics ' ,  
 C    &Nr, 1, myThid )  
 C     CALL PLOT_FIELD_XYZRL( gS, ' GS exiting dyanmics ' ,  
 C    &Nr, 1, myThid )  
 C     CALL PLOT_FIELD_XYZRL( gT, ' GT exiting dyanmics ' ,  
 C    &Nr, 1, myThid )  
 C     CALL PLOT_FIELD_XYZRL( phiHyd, ' phiHyd exiting dyanmics ' ,  
 C    &Nr, 1, myThid )  
   
   
682        RETURN        RETURN
683        END        END
Legend:
Removed from v.1.49  
changed lines
  Added in v.1.65
  ViewVC Help
Powered by ViewVC 1.1.22