/[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.55 by heimbach, Mon Jan 8 16:37:43 2001 UTC
# Line 20  C     | ===== Line 20  C     | =====
20  C     | C*P* comments indicating place holders for which code is |  C     | C*P* comments indicating place holders for which code is |
21  C     |      presently being developed.                          |  C     |      presently being developed.                          |
22  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.  
   
23        IMPLICIT NONE        IMPLICIT NONE
24    
25  C     == Global variables ===  C     == Global variables ===
# Line 36  C     == Global variables === Line 30  C     == Global variables ===
30  #include "DYNVARS.h"  #include "DYNVARS.h"
31  #include "GRID.h"  #include "GRID.h"
32    
 #ifdef ALLOW_KPP  
 #include "KPPMIX.h"  
 #endif  
   
33  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
34  #include "tamc.h"  # include "tamc.h"
35  #include "tamc_keys.h"  # include "tamc_keys.h"
36    #endif /* ALLOW_AUTODIFF_TAMC */
37    
38    #ifdef ALLOW_KPP
39    # include "KPP.h"
40  #endif  #endif
41    
42  C     == Routine arguments ==  C     == Routine arguments ==
# Line 89  C                      surface height Line 83  C                      surface height
83  C                      anomaly.  C                      anomaly.
84  C     etaSurfX,      - Holds surface elevation gradient in X and Y.  C     etaSurfX,      - Holds surface elevation gradient in X and Y.
85  C     etaSurfY  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.  
86  C     KappaRT,       - Total diffusion in vertical for T and S.  C     KappaRT,       - Total diffusion in vertical for T and S.
87  C     KappaRS          (background + spatially varying, isopycnal term).  C     KappaRS          (background + spatially varying, isopycnal term).
88  C     iMin, iMax     - Ranges and sub-block indices on which calculations  C     iMin, iMax     - Ranges and sub-block indices on which calculations
89  C     jMin, jMax       are applied.  C     jMin, jMax       are applied.
90  C     bi, bj  C     bi, bj
91  C     k, kUp,        - Index for layer above and below. kUp and kDown  C     k, kup,        - Index for layer above and below. kup and kDown
92  C     kDown, kM1       are switched with layer to be the appropriate  C     kDown, km1       are switched with layer to be the appropriate
93  C                      index into fVerTerm.  C                      index into fVerTerm.
94        _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
95        _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 128  C                      index into fVerTe Line 119  C                      index into fVerTe
119        _RL rhotmp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhotmp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
120        _RL etaSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL etaSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
121        _RL etaSurfY(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)  
122        _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
123        _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
124        _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
125        _RL KappaRV (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRV (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
126          _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
127          _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
128          _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
129    
130  #ifdef INCLUDE_CONVECT_CALL  C This is currently also used by IVDC and Diagnostics
131    C #ifdef INCLUDE_CONVECT_CALL
132        _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
133  #endif  C #endif
134    
135        INTEGER iMin, iMax        INTEGER iMin, iMax
136        INTEGER jMin, jMax        INTEGER jMin, jMax
137        INTEGER bi, bj        INTEGER bi, bj
138        INTEGER i, j        INTEGER i, j
139        INTEGER k, kM1, kUp, kDown        INTEGER k, km1, kup, kDown
140        LOGICAL BOTTOM_LAYER        LOGICAL BOTTOM_LAYER
141    
142  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 154  C                      index into fVerTe Line 145  C                      index into fVerTe
145    
146        INTEGER act1, act2, act3, act4        INTEGER act1, act2, act3, act4
147        INTEGER max1, max2, max3        INTEGER max1, max2, max3
148        INTEGER ikact, iikey,kkey        INTEGER iikey, kkey
149        INTEGER maximpl        INTEGER maximpl
150  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
151    
152  C---    The algorithm...  C---    The algorithm...
153  C  C
# Line 171  C Line 162  C
162  C       "Calculation of Gs"  C       "Calculation of Gs"
163  C       ===================  C       ===================
164  C       This is where all the accelerations and tendencies (ie.  C       This is where all the accelerations and tendencies (ie.
165  C       phiHydysics, parameterizations etc...) are calculated  C       physics, parameterizations etc...) are calculated
166  C         rVel = sum_r ( div. u[n] )  C         rVel = sum_r ( div. u[n] )
167  C         rho = rho ( theta[n], salt[n] )  C         rho = rho ( theta[n], salt[n] )
168  C         b   = b(rho, theta)  C         b   = b(rho, theta)
# Line 206  C--- Line 197  C---
197  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
198  C--   dummy statement to end declaration part  C--   dummy statement to end declaration part
199        ikey = 1        ikey = 1
200  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
201    
202    
203  C--   Set up work arrays with valid (i.e. not NaN) values  C--   Set up work arrays with valid (i.e. not NaN) values
204  C     These inital values do not alter the numerical results. They  C     These inital values do not alter the numerical results. They
# Line 226  C     uninitialised but inert locations. Line 218  C     uninitialised but inert locations.
218          pTerm(i,j)   = 0. _d 0          pTerm(i,j)   = 0. _d 0
219          fZon(i,j)    = 0. _d 0          fZon(i,j)    = 0. _d 0
220          fMer(i,j)    = 0. _d 0          fMer(i,j)    = 0. _d 0
221          DO K=1,Nr          DO k=1,Nr
222           phiHyd (i,j,k)  = 0. _d 0           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  
223           KappaRU(i,j,k) = 0. _d 0           KappaRU(i,j,k) = 0. _d 0
224           KappaRV(i,j,k) = 0. _d 0           KappaRV(i,j,k) = 0. _d 0
225             sigmaX(i,j,k) = 0. _d 0
226             sigmaY(i,j,k) = 0. _d 0
227             sigmaR(i,j,k) = 0. _d 0
228          ENDDO          ENDDO
229          rhoKM1 (i,j) = 0. _d 0          rhoKM1 (i,j) = 0. _d 0
230          rhok   (i,j) = 0. _d 0          rhok   (i,j) = 0. _d 0
# Line 247  C     uninitialised but inert locations. Line 239  C     uninitialised but inert locations.
239    
240  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
241  C--   HPF directive to help TAMC  C--   HPF directive to help TAMC
242  !HPF$ INDEPENDENT  CHPF$ INDEPENDENT
243  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
244    
245        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
246    
247  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
248  C--    HPF directive to help TAMC  C--    HPF directive to help TAMC
249  !HPF$  INDEPENDENT, NEW (rTrans,rVel,fVerT,fVerS,fVerU,fVerV  CHPF$  INDEPENDENT, NEW (rTrans,rVel,fVerT,fVerS,fVerU,fVerV
250  !HPF$&                  ,phiHyd,K13,K23,K33,KapGM  CHPF$&                  ,phiHyd,utrans,vtrans,maskc,xA,yA
251  !HPF$&                  ,utrans,vtrans,maskc,xA,yA  CHPF$&                  ,KappaRT,KappaRS,KappaRU,KappaRV
252  !HPF$&                  ,KappaRT,KappaRS,KappaRU,KappaRV  CHPF$&                  )
253  !HPF$&                  )  #endif /* ALLOW_AUTODIFF_TAMC */
 #endif  
254    
255         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
256    
# Line 278  C--    HPF directive to help TAMC Line 269  C--    HPF directive to help TAMC
269            ikey = (act1 + 1) + act2*max1            ikey = (act1 + 1) + act2*max1
270       &                      + act3*max1*max2       &                      + act3*max1*max2
271       &                      + act4*max1*max2*max3       &                      + act4*max1*max2*max3
272  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
273    
274  C--     Set up work arrays that need valid initial values  C--     Set up work arrays that need valid initial values
275          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
# Line 295  C--     Set up work arrays that need val Line 286  C--     Set up work arrays that need val
286            fVerV (i,j,1) = 0. _d 0            fVerV (i,j,1) = 0. _d 0
287            fVerV (i,j,2) = 0. _d 0            fVerV (i,j,2) = 0. _d 0
288            phiHyd(i,j,1) = 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  
289           ENDDO           ENDDO
290          ENDDO          ENDDO
291    
# Line 319  C--     Set up work arrays that need val Line 306  C--     Set up work arrays that need val
306          jMin = 1-OLy+1          jMin = 1-OLy+1
307          jMax = sNy+OLy          jMax = sNy+OLy
308    
309            k = 1
310          K = 1          BOTTOM_LAYER = k .EQ. Nr
         BOTTOM_LAYER = K .EQ. Nr  
311    
312  #ifdef DO_PIPELINED_CORRECTION_STEP  #ifdef DO_PIPELINED_CORRECTION_STEP
313  C--     Calculate gradient of surface pressure  C--     Calculate gradient of surface pressure
# Line 331  C--     Calculate gradient of surface pr Line 317  C--     Calculate gradient of surface pr
317       I       myThid)       I       myThid)
318  C--     Update fields in top level according to tendency terms  C--     Update fields in top level according to tendency terms
319          CALL CORRECTION_STEP(          CALL CORRECTION_STEP(
320       I       bi,bj,iMin,iMax,jMin,jMax,K,       I       bi,bj,iMin,iMax,jMin,jMax,k,
321       I       etaSurfX,etaSurfY,myTime,myThid)       I       etaSurfX,etaSurfY,myTime,myThid)
322    
323  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
324          IF (openBoundaries) THEN          IF (openBoundaries) THEN
325  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
326  CADJ STORE uvel (:,:,k,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  CADJ STORE uvel (:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
327  CADJ STORE vvel (:,:,k,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  CADJ STORE vvel (:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
328  CADJ STORE theta(:,:,k,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
329  CADJ STORE salt(:,:,k,bi,bj)   = comlev1_2d, key = ikey, byte = isbyte  CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
330  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
331             CALL APPLY_OBCS1( bi, bj, K, myThid )             CALL APPLY_OBCS1( bi, bj, k, myThid )
332          END IF          END IF
333  #endif  #endif
334    
335          IF ( .NOT. BOTTOM_LAYER ) THEN          IF ( .NOT. BOTTOM_LAYER ) THEN
336  C--      Update fields in layer below according to tendency terms  C--      Update fields in layer below according to tendency terms
337           CALL CORRECTION_STEP(           CALL CORRECTION_STEP(
338       I        bi,bj,iMin,iMax,jMin,jMax,K+1,       I        bi,bj,iMin,iMax,jMin,jMax,k+1,
339       I        etaSurfX,etaSurfY,myTime,myThid)       I        etaSurfX,etaSurfY,myTime,myThid)
340  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
341           IF (openBoundaries) THEN           IF (openBoundaries) THEN
342  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
343  CADJ STORE uvel (:,:,k,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  CADJ STORE uvel (:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
344  CADJ STORE vvel (:,:,k,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  CADJ STORE vvel (:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
345  CADJ STORE theta(:,:,k,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
346  CADJ STORE salt(:,:,k,bi,bj)   = comlev1_2d, key = ikey, byte = isbyte  CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
347  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
348              CALL APPLY_OBCS1( bi, bj, K+1, myThid )              CALL APPLY_OBCS1( bi, bj, k+1, myThid )
349           END IF           END IF
350  #endif  #endif
351          ENDIF          ENDIF
352  #endif  #endif
353    
354  C--     Density of 1st level (below W(1)) reference to level 1  C--     Density of 1st level (below W(1)) reference to level 1
355  #ifdef  INCLUDE_FIND_RHO_CALL  #ifdef  INCLUDE_FIND_RHO_CALL
356  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
357  CADJ STORE theta(:,:,k,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
358  CADJ STORE salt (:,:,k,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
359  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
360          CALL FIND_RHO(          CALL FIND_RHO(
361       I     bi, bj, iMin, iMax, jMin, jMax, K, K, eosType,       I     bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,
362       O     rhoKm1,       O     rhoKm1,
363       I     myThid )       I     myThid )
364  #endif  #endif
365    
366          IF (       (.NOT. BOTTOM_LAYER)          IF (.NOT. BOTTOM_LAYER) THEN
367  #ifdef ALLOW_KPP  
      &       .AND. (.NOT.usingKPPmixing) ! CONVECT not needed with KPP mixing  
 #endif  
      &     ) THEN  
368  C--      Check static stability with layer below  C--      Check static stability with layer below
369  C--      and mix as needed.  C--      and mix as needed.
370  #ifdef  INCLUDE_FIND_RHO_CALL  #ifdef  INCLUDE_FIND_RHO_CALL
371  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
372  CADJ STORE theta(:,:,k,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  CADJ STORE theta(:,:,k+1,bi,bj) = comlev1_bibj
373  CADJ STORE salt (:,:,k,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  CADJ &   , key = ikey, byte = isbyte
374  #endif  CADJ STORE salt (:,:,k+1,bi,bj) = comlev1_bibj
375    CADJ &   , key = ikey, byte = isbyte
376    #endif /* ALLOW_AUTODIFF_TAMC */
377           CALL FIND_RHO(           CALL FIND_RHO(
378       I      bi, bj, iMin, iMax, jMin, jMax, K+1, K, eosType,       I      bi, bj, iMin, iMax, jMin, jMax, k+1, k, eosType,
379       O      rhoKp1,       O      rhoKp1,
380       I      myThid )       I      myThid )
381  #endif  #endif
382    
383    #ifdef ALLOW_AUTODIFF_TAMC
384    CADJ STORE rhoKm1(:,:)  = comlev1_bibj, key = ikey, byte = isbyte
385    CADJ STORE rhoKp1(:,:)  = comlev1_bibj, key = ikey, byte = isbyte
386    #endif /* ALLOW_AUTODIFF_TAMC */
387    
388  #ifdef  INCLUDE_CONVECT_CALL  #ifdef  INCLUDE_CONVECT_CALL
389    
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE rhoKm1(:,:)  = comlev1_2d, key = ikey, byte = isbyte  
 CADJ STORE rhoKp1(:,:)  = comlev1_2d, key = ikey, byte = isbyte  
 #endif  
390           CALL CONVECT(           CALL CONVECT(
391       I       bi,bj,iMin,iMax,jMin,jMax,K+1,rhoKm1,rhoKp1,       I       bi,bj,iMin,iMax,jMin,jMax,k+1,rhoKm1,rhoKp1,
392       U       ConvectCount,       U       ConvectCount,
393       I       myTime,myIter,myThid)       I       myTime,myIter,myThid)
394    
395  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
396  CADJ STORE theta(:,:,k+1,bi,bj),theta(:,:,k,bi,bj)  CADJ STORE theta(:,:,k+1,bi,bj),theta(:,:,k,bi,bj)
397  CADJ &     = comlev1_2d, key = ikey, byte = isbyte  CADJ &     = comlev1_bibj, key = ikey, byte = isbyte
398  CADJ STORE salt (:,:,k+1,bi,bj),salt (:,:,k,bi,bj)  CADJ STORE salt (:,:,k+1,bi,bj),salt (:,:,k,bi,bj)
399  CADJ &     = comlev1_2d, key = ikey, byte = isbyte  CADJ &     = comlev1_bibj, key = ikey, byte = isbyte
400  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
401    
402  #endif  #endif
403    
404  C--      Implicit Vertical Diffusion for Convection  C--      Implicit Vertical Diffusion for Convection
405           IF (ivdc_kappa.NE.0.) CALL CALC_IVDC(           IF (ivdc_kappa.NE.0.) THEN
406       I       bi,bj,iMin,iMax,jMin,jMax,K+1,rhoKm1,rhoKp1,              CALL CALC_IVDC(
407         I       bi,bj,iMin,iMax,jMin,jMax,k+1,rhoKm1,rhoKp1,
408       U       ConvectCount, KappaRT, KappaRS,       U       ConvectCount, KappaRT, KappaRS,
409       I       myTime,myIter,myThid)       I       myTime,myIter,myThid)
410  CRG: do we need do store STORE KappaRT, KappaRS ?           ENDIF
411    
412  C--      Recompute density after mixing  C--      Recompute density after mixing
413  #ifdef  INCLUDE_FIND_RHO_CALL  #ifdef  INCLUDE_FIND_RHO_CALL
414           CALL FIND_RHO(           CALL FIND_RHO(
415       I      bi, bj, iMin, iMax, jMin, jMax, K, K, eosType,       I      bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,
416       O      rhoKm1,       O      rhoKm1,
417       I      myThid )       I      myThid )
418  #endif  #endif
419          ENDIF          ENDIF
420    
421  C--     Calculate buoyancy  C--     Calculate buoyancy
422          CALL CALC_BUOYANCY(          CALL CALC_BUOYANCY(
423       I      bi,bj,iMin,iMax,jMin,jMax,K,rhoKm1,       I      bi,bj,iMin,iMax,jMin,jMax,k,rhoKm1,
424       O      buoyKm1,       O      buoyKm1,
425       I      myThid )       I      myThid )
426    
427  C--     Integrate hydrostatic balance for phiHyd with BC of  C--     Integrate hydrostatic balance for phiHyd with BC of
428  C--     phiHyd(z=0)=0  C--     phiHyd(z=0)=0
429          CALL CALC_PHI_HYD(          CALL CALC_PHI_HYD(
430       I      bi,bj,iMin,iMax,jMin,jMax,K,buoyKm1,buoyKm1,       I      bi,bj,iMin,iMax,jMin,jMax,k,buoyKm1,buoyKm1,
431       U      phiHyd,       U      phiHyd,
432       I      myThid )       I      myThid )
433    
434  C----------------------------------------------  #ifdef ALLOW_GMREDI
435  C--     start of downward loop          IF ( useGMRedi ) THEN
436  C----------------------------------------------          CALL GRAD_SIGMA(
437          DO K=2,Nr       I            bi, bj, iMin, iMax, jMin, jMax, k,
438         I            rhoKm1, rhoKm1, rhoKm1,
439         O            sigmaX, sigmaY, sigmaR,
440         I            myThid )
441            ELSE
442             DO j=1-OLy,sNy+OLy
443              DO i=1-OLx,sNx+OLx
444               sigmaX(i,j,k) = 0. _d 0
445               sigmaY(i,j,k) = 0. _d 0
446               sigmaR(i,j,k) = 0. _d 0
447              ENDDO
448             ENDDO
449            ENDIF
450    #endif
451    
452    C--     Start of downward loop
453            DO k=2,Nr
454    
455  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
456           kkey = ikact*(Nr-2+1) + (k-2) + 1           kkey = (ikey-1)*(Nr-2+1) + (k-2) + 1
457  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
458    
459           BOTTOM_LAYER = K .EQ. Nr           BOTTOM_LAYER = k .EQ. Nr
460    
461  #ifdef DO_PIPELINED_CORRECTION_STEP  #ifdef DO_PIPELINED_CORRECTION_STEP
462           IF ( .NOT. BOTTOM_LAYER ) THEN           IF ( .NOT. BOTTOM_LAYER ) THEN
463  C--       Update fields in layer below according to tendency terms  C--       Update fields in layer below according to tendency terms
464            CALL CORRECTION_STEP(            CALL CORRECTION_STEP(
465       I         bi,bj,iMin,iMax,jMin,jMax,K+1,       I         bi,bj,iMin,iMax,jMin,jMax,k+1,
466       I         etaSurfX,etaSurfY,myTime,myThid)       I         etaSurfX,etaSurfY,myTime,myThid)
467  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
468            IF (openBoundaries) THEN            IF (openBoundaries) THEN
469  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
470  CADJ STORE uvel (:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  CADJ STORE uvel (:,:,k,bi,bj) = comlev1_bibj_k
471  CADJ STORE vvel (:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  CADJ &   , key = kkey, byte = isbyte
472  CADJ STORE theta(:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  CADJ STORE vvel (:,:,k,bi,bj) = comlev1_bibj_k
473  CADJ STORE salt(:,:,k,bi,bj)   = comlev1_2d, key = ikey, byte = isbyte  CADJ &   , key = kkey, byte = isbyte
474  #endif  CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k
475               CALL APPLY_OBCS1( bi, bj, K+1, myThid )  CADJ &   , key = kkey, byte = isbyte
476    CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k
477    CADJ &   , key = kkey, byte = isbyte
478    #endif /* ALLOW_AUTODIFF_TAMC */
479                 CALL APPLY_OBCS1( bi, bj, k+1, myThid )
480            END IF            END IF
481  #endif  #endif
482           ENDIF           ENDIF
483  #endif  #endif /* DO_PIPELINED_CORRECTION_STEP */
484    
485  C--      Density of K level (below W(K)) reference to K level  C--      Density of k level (below W(k)) reference to k level
486  #ifdef  INCLUDE_FIND_RHO_CALL  #ifdef  INCLUDE_FIND_RHO_CALL
487  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
488  CADJ STORE theta(:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k
489  CADJ STORE salt (:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  CADJ &   , key = kkey, byte = isbyte
490  #endif  CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k
491    CADJ &   , key = kkey, byte = isbyte
492    #endif /* ALLOW_AUTODIFF_TAMC */
493           CALL FIND_RHO(           CALL FIND_RHO(
494       I      bi, bj, iMin, iMax, jMin, jMax,  K, K, eosType,       I      bi, bj, iMin, iMax, jMin, jMax,  k, k, eosType,
495       O      rhoK,       O      rhoK,
496       I      myThid )       I      myThid )
497    
498    #ifdef ALLOW_AUTODIFF_TAMC
499    cph(   storing not necessary
500    cphCADJ STORE rhoK(:,:)  = comlev1_bibj_k, key = kkey, byte = isbyte
501    cph)
502    #endif /* ALLOW_AUTODIFF_TAMC */
503  #endif  #endif
504           IF (       (.NOT. BOTTOM_LAYER)  
505  #ifdef ALLOW_KPP           IF (.NOT. BOTTOM_LAYER) THEN
506       &       .AND. (.NOT.usingKPPmixing) ! CONVECT not needed with KPP mixing  
 #endif  
      &      ) THEN  
507  C--       Check static stability with layer below and mix as needed.  C--       Check static stability with layer below and mix as needed.
508  C--       Density of K+1 level (below W(K+1)) reference to K level.  C--       Density of k+1 level (below W(k+1)) reference to k level.
509  #ifdef  INCLUDE_FIND_RHO_CALL  #ifdef  INCLUDE_FIND_RHO_CALL
510  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
511  CADJ STORE theta(:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  CADJ STORE theta(:,:,k+1,bi,bj) = comlev1_bibj_k
512  CADJ STORE salt (:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  CADJ &   , key = kkey, byte = isbyte
513  #endif  CADJ STORE salt (:,:,k+1,bi,bj) = comlev1_bibj_k
514    CADJ &   , key = kkey, byte = isbyte
515    #endif /* ALLOW_AUTODIFF_TAMC */
516            CALL FIND_RHO(            CALL FIND_RHO(
517       I       bi, bj, iMin, iMax, jMin, jMax,  K+1, K, eosType,       I       bi, bj, iMin, iMax, jMin, jMax,  k+1, k, eosType,
518       O       rhoKp1,       O       rhoKp1,
519       I       myThid )       I       myThid )
 #endif  
   
520  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
521  CADJ STORE rhok  (:,:)   = comlev1_3d, key = kkey, byte = isbyte  CADJ STORE rhoKp1(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
522  CADJ STORE rhoKm1(:,:)   = comlev1_3d, key = kkey, byte = isbyte  #endif /* ALLOW_AUTODIFF_TAMC */
 CADJ STORE rhoKp1(:,:)   = comlev1_3d, key = kkey, byte = isbyte  
523  #endif  #endif
524    
525  #ifdef  INCLUDE_CONVECT_CALL  #ifdef  INCLUDE_CONVECT_CALL
526            CALL CONVECT(            CALL CONVECT(
527       I        bi,bj,iMin,iMax,jMin,jMax,K+1,rhoK,rhoKp1,       I        bi,bj,iMin,iMax,jMin,jMax,k+1,rhoK,rhoKp1,
528       U        ConvectCount,       U        ConvectCount,
529       I        myTime,myIter,myThid)       I        myTime,myIter,myThid)
530  #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE theta(:,:,k+1,bi,bj),theta(:,:,k,bi,bj)  
 CADJ &     = comlev1_3d, key = kkey, byte = isbyte  
 CADJ STORE salt (:,:,k+1,bi,bj),salt (:,:,k,bi,bj)  
 CADJ &     = comlev1_3d, key = kkey, byte = isbyte  
 #endif  
531  #endif  #endif
532    
533  C--      Implicit Vertical Diffusion for Convection  C--      Implicit Vertical Diffusion for Convection
534           IF (ivdc_kappa.NE.0.) THEN           IF (ivdc_kappa.NE.0.) THEN
535    #ifdef ALLOW_AUTODIFF_TAMC
536    CADJ STORE rhoKm1(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
537    #endif /* ALLOW_AUTODIFF_TAMC */
538              CALL CALC_IVDC(              CALL CALC_IVDC(
539       I       bi,bj,iMin,iMax,jMin,jMax,K+1,rhoKm1,rhoKp1,       I       bi,bj,iMin,iMax,jMin,jMax,k+1,rhoKm1,rhoKp1,
540       U       ConvectCount, KappaRT, KappaRS,       U       ConvectCount, KappaRT, KappaRS,
541       I       myTime,myIter,myThid)       I       myTime,myIter,myThid)
 CRG: do we need do store STORE KappaRT, KappaRS ?  
542           END IF           END IF
543    
544  C--       Recompute density after mixing  C--       Recompute density after mixing
545  #ifdef  INCLUDE_FIND_RHO_CALL  #ifdef  INCLUDE_FIND_RHO_CALL
546    #ifdef ALLOW_AUTODIFF_TAMC
547    CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k
548    CADJ &   , key = kkey, byte = isbyte
549    CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k
550    CADJ &   , key = kkey, byte = isbyte
551    #endif /* ALLOW_AUTODIFF_TAMC */
552            CALL FIND_RHO(            CALL FIND_RHO(
553       I       bi, bj, iMin, iMax, jMin, jMax, K, K, eosType,       I       bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,
554       O       rhoK,       O       rhoK,
555       I       myThid )       I       myThid )
556  #endif  #endif
557    
558    C--            IF (.NOT. BOTTOM_LAYER) ends here
559           ENDIF           ENDIF
560    
561  C--      Calculate buoyancy  C--      Calculate buoyancy
562           CALL CALC_BUOYANCY(           CALL CALC_BUOYANCY(
563       I       bi,bj,iMin,iMax,jMin,jMax,K,rhoK,       I       bi,bj,iMin,iMax,jMin,jMax,k,rhoK,
564       O       buoyK,       O       buoyK,
565       I       myThid )       I       myThid )
566    
567  C--      Integrate hydrostatic balance for phiHyd with BC of  C--      Integrate hydrostatic balance for phiHyd with BC of
568  C--      phiHyd(z=0)=0  C--      phiHyd(z=0)=0
569           CALL CALC_PHI_HYD(           CALL CALC_PHI_HYD(
570       I        bi,bj,iMin,iMax,jMin,jMax,K,buoyKm1,buoyK,       I        bi,bj,iMin,iMax,jMin,jMax,k,buoyKm1,buoyK,
571       U        phiHyd,       U        phiHyd,
572       I        myThid )       I        myThid )
573  C--      Calculate iso-neutral slopes for the GM/Redi parameterisation  
574  #ifdef  INCLUDE_FIND_RHO_CALL  #ifdef  INCLUDE_FIND_RHO_CALL
575    C--      Calculate iso-neutral slopes for the GM/Redi parameterisation
576    
577    #ifdef ALLOW_AUTODIFF_TAMC
578    CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k
579    CADJ &   , key = kkey, byte = isbyte
580    CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k
581    CADJ &   , key = kkey, byte = isbyte
582    #endif /* ALLOW_AUTODIFF_TAMC */
583    
584           CALL FIND_RHO(           CALL FIND_RHO(
585       I        bi, bj, iMin, iMax, jMin, jMax, K-1, K, eosType,       I        bi, bj, iMin, iMax, jMin, jMax, k-1, k, eosType,
586       O        rhoTmp,       O        rhoTmp,
587       I        myThid )       I        myThid )
588  #endif  #endif
589    
590  #ifdef  INCLUDE_CALC_ISOSLOPES_CALL  
591  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_GMREDI
592  CADJ STORE rhoTmp(:,:)  = comlev1_3d, key = kkey, byte = isbyte           IF ( useGMRedi ) THEN
593  CADJ STORE rhok  (:,:)  = comlev1_3d, key = kkey, byte = isbyte           CALL GRAD_SIGMA(
594  CADJ STORE rhoKm1(:,:)  = comlev1_3d, key = kkey, byte = isbyte       I             bi, bj, iMin, iMax, jMin, jMax, k,
595  CADJ STORE kapgm (:,:)  = comlev1_3d, key = kkey, byte = isbyte       I             rhoK, rhotmp, rhoK,
596  #endif       O             sigmaX, sigmaY, sigmaR,
597           CALL CALC_ISOSLOPES(       I             myThid )
598       I        bi, bj, iMin, iMax, jMin, jMax, K,           ELSE
599       I        rhoKm1, rhoK, rhotmp,            DO j=1-OLy,sNy+OLy
600       O        K13, K23, K33, KapGM,             DO i=1-OLx,sNx+OLx
601       I        myThid )              sigmaX(i,j,k) = 0. _d 0
602                sigmaY(i,j,k) = 0. _d 0
603                sigmaR(i,j,k) = 0. _d 0
604               ENDDO
605              ENDDO
606             ENDIF
607  #endif  #endif
608    
609           DO J=jMin,jMax           DO J=jMin,jMax
# Line 577  CADJ STORE kapgm (:,:)  = comlev1_3d, ke Line 614  CADJ STORE kapgm (:,:)  = comlev1_3d, ke
614             buoyKm1(I,J) = buoyK(I,J)             buoyKm1(I,J) = buoyK(I,J)
615            ENDDO            ENDDO
616           ENDDO           ENDDO
617          ENDDO  
618  C--     end of k loop  C--     end of k loop
619            ENDDO
620    
621    C     Determines forcing terms based on external fields
622    C     relaxation terms, etc.
623          CALL EXTERNAL_FORCING_SURF(
624         I             bi, bj, iMin, iMax, jMin, jMax,
625         I             myThid )
626    
627    #ifdef ALLOW_AUTODIFF_TAMC
628    
629    CADJ STORE surfacetendencyu(:,:,bi,bj)
630    CADJ &   , surfacetendencyv(:,:,bi,bj)
631    CADJ &   , surfacetendencys(:,:,bi,bj)
632    CADJ &   , surfacetendencyt(:,:,bi,bj)
633    CADJ &                        = comlev1_bibj, key=ikey, byte=isbyte
634    
635    # ifdef ALLOW_GMREDI
636    CADJ STORE sigmaX(:,:,:) = comlev1, key=ikey, byte=isbyte
637    CADJ STORE sigmaY(:,:,:) = comlev1, key=ikey, byte=isbyte
638    CADJ STORE sigmaR(:,:,:) = comlev1, key=ikey, byte=isbyte
639    # endif /* ALLOW_GMREDI */
640    
641    #endif /* ALLOW_AUTODIFF_TAMC */
642    
643    #ifdef ALLOW_GMREDI
644            IF (useGMRedi) THEN
645              DO k=1, Nr
646                CALL GMREDI_CALC_TENSOR(
647         I             bi, bj, iMin, iMax, jMin, jMax, k,
648         I             sigmaX, sigmaY, sigmaR,
649         I             myThid )
650              ENDDO
651  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
652  CADJ STORE theta(:,:,:,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte          ELSE
653  CADJ STORE salt (:,:,:,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte            DO k=1, Nr
654  CADJ STORE uvel (:,:,:,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte              CALL GMREDI_CALC_TENSOR_DUMMY(
655  CADJ STORE vvel (:,:,:,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte       I             bi, bj, iMin, iMax, jMin, jMax, k,
656         I             sigmaX, sigmaY, sigmaR,
657         I             myThid )
658              ENDDO
659    #endif /* ALLOW_AUTODIFF_TAMC */
660            ENDIF
661  #endif  #endif
662    
663    #ifdef ALLOW_AUTODIFF_TAMC
664    CADJ STORE KappaRT(:,:,:)     = comlev1_bibj, key=ikey, byte=isbyte
665    CADJ STORE KappaRS(:,:,:)     = comlev1_bibj, key=ikey, byte=isbyte
666    
667    C-- R.G. We need to define a new tape since Kw use mythid instead of bi,bj
668    CADJ STORE Kwx(:,:,:,myThid)  = comlev1_bibj, key=ikey, byte=isbyte
669    CADJ STORE Kwy(:,:,:,myThid)  = comlev1_bibj, key=ikey, byte=isbyte
670    CADJ STORE Kwz(:,:,:,myThid)  = comlev1_bibj, key=ikey, byte=isbyte
671    
672    CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
673    CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
674    CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
675    CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
676    
677    C--     dummy initialization to break data flow because
678    C--     calc_div_ghat has a condition for initialization
679            DO J=jMin,jMax
680               DO I=iMin,iMax
681                  cg2d_b(i,j,bi,bj) = 0.0
682               ENDDO
683            ENDDO
684    #endif /* ALLOW_AUTODIFF_TAMC */
685    
686  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
687  C----------------------------------------------  C--   Compute KPP mixing coefficients
688  C--     Compute KPP mixing coefficients          IF (useKPP) THEN
689  C----------------------------------------------  
690          IF (usingKPPmixing) THEN            CALL TIMER_START('KPP_CALC               [DYNAMICS]', myThid)
691  #ifdef ALLOW_AUTODIFF_TAMC            CALL KPP_CALC(
692  CADJ STORE fu  (:,:  ,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte       I                  bi, bj, myTime, myThid )
693  CADJ STORE fv  (:,:  ,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte            CALL TIMER_STOP ('KPP_CALC               [DYNAMICS]', myThid)
694  #endif  
695           CALL TIMER_START('KVMIX (FIND KPP COEFFICIENTS) [DYNAMICS]'  #ifdef ALLOW_AUTODIFF_TAMC
696       I          , myThid)          ELSE
697           CALL KVMIX(            DO j=1-OLy,sNy+OLy
698       I               bi, bj, myTime, myThid )              DO i=1-OLx,sNx+OLx
699           CALL TIMER_STOP ('KVMIX (FIND KPP COEFFICIENTS) [DYNAMICS]'                KPPhbl (i,j,bi,bj) = 1.0
700       I        , myThid)                KPPfrac(i,j,bi,bj) = 0.0
701                  DO k = 1,Nr
702                     KPPghat   (i,j,k,bi,bj) = 0.0
703                     KPPviscAz (i,j,k,bi,bj) = viscAz
704                     KPPdiffKzT(i,j,k,bi,bj) = diffKzT
705                     KPPdiffKzS(i,j,k,bi,bj) = diffKzS
706                  ENDDO
707                ENDDO
708              ENDDO
709    #endif /* ALLOW_AUTODIFF_TAMC */
710          ENDIF          ENDIF
 #endif  
711    
712  C----------------------------------------------  #ifdef ALLOW_AUTODIFF_TAMC
713  C--     start of upward loop  CADJ STORE KPPghat   (:,:,:,bi,bj)
714  C----------------------------------------------  CADJ &   , KPPviscAz (:,:,:,bi,bj)
715          DO K = Nr, 1, -1  CADJ &   , KPPdiffKzT(:,:,:,bi,bj)
716    CADJ &   , KPPdiffKzS(:,:,:,bi,bj)
717           kM1  =max(1,k-1)   ! Points to level above k (=k-1)  CADJ &   , KPPfrac   (:,:  ,bi,bj)
718           kUp  =1+MOD(k+1,2) ! Cycles through 1,2 to point to layer above  CADJ &                 = comlev1_bibj, key=ikey, byte=isbyte
719           kDown=1+MOD(k,2)   ! Cycles through 2,1 to point to current layer  #endif /* ALLOW_AUTODIFF_TAMC */
720    
721    #endif /* ALLOW_KPP */
722    
723    C--     Start of upward loop
724            DO k = Nr, 1, -1
725    
726    C--      km1    Points to level above k (=k-1)
727    C--      kup    Cycles through 1,2 to point to layer above
728    C--      kDown  Cycles through 2,1 to point to current layer
729    
730             km1  =max(1,k-1)
731             kup  =1+MOD(k+1,2)
732             kDown=1+MOD(k,2)
733    
734           iMin = 1-OLx+2           iMin = 1-OLx+2
735           iMax = sNx+OLx-1           iMax = sNx+OLx-1
# Line 620  C--------------------------------------- Line 737  C---------------------------------------
737           jMax = sNy+OLy-1           jMax = sNy+OLy-1
738    
739  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
740           kkey = ikact*(Nr-1+1) + (k-1) + 1           kkey = (ikey-1)*(Nr-1+1) + (k-1) + 1
 #endif  
741    
742  #ifdef ALLOW_AUTODIFF_TAMC  CADJ STORE rvel  (:,:,kDown) = comlev1_bibj_k, key = kkey, byte = isbyte
743  CADJ STORE rvel  (:,:,kDown)  = comlev1_3d, key = kkey, byte = isbyte  CADJ STORE rTrans(:,:)       = comlev1_bibj_k, key = kkey, byte = isbyte
744  CADJ STORE rTrans(:,:)        = comlev1_3d, key = kkey, byte = isbyte  CADJ STORE KappaRT(:,:,:)    = comlev1_bibj_k, key = kkey, byte = isbyte
745  CADJ STORE KappaRT(:,:,:)     = comlev1_3d, key = kkey, byte = isbyte  CADJ STORE KappaRS(:,:,:)    = comlev1_bibj_k, key = kkey, byte = isbyte
746  CADJ STORE KappaRS(:,:,:)     = comlev1_3d, key = kkey, byte = isbyte  #endif /* ALLOW_AUTODIFF_TAMC */
 #endif  
747    
748  C--      Get temporary terms used by tendency routines  C--      Get temporary terms used by tendency routines
749           CALL CALC_COMMON_FACTORS (           CALL CALC_COMMON_FACTORS (
750       I        bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,       I        bi,bj,iMin,iMax,jMin,jMax,k,km1,kup,kDown,
751       O        xA,yA,uTrans,vTrans,rTrans,rVel,maskC,maskUp,       O        xA,yA,uTrans,vTrans,rTrans,rVel,maskC,maskUp,
752       I        myThid)       I        myThid)
753    
754  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
755          IF (openBoundaries) THEN          IF (openBoundaries) THEN
756           CALL APPLY_OBCS3( bi, bj, K, Kup, rTrans, rVel, myThid )           CALL APPLY_OBCS3( bi, bj, k, kup, rTrans, rVel, myThid )
757          ENDIF          ENDIF
758  #endif  #endif
759    
760  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL
761  C--      Calculate the total vertical diffusivity  C--      Calculate the total vertical diffusivity
762           CALL CALC_DIFFUSIVITY(           CALL CALC_DIFFUSIVITY(
763       I        bi,bj,iMin,iMax,jMin,jMax,K,       I        bi,bj,iMin,iMax,jMin,jMax,k,
764       I        maskC,maskUp,KapGM,K33,       I        maskC,maskUp,
765       O        KappaRT,KappaRS,KappaRU,KappaRV,       O        KappaRT,KappaRS,KappaRU,KappaRV,
766       I        myThid)       I        myThid)
767  #endif  #endif
768  C--      Calculate accelerations in the momentum equations  C--      Calculate accelerations in the momentum equations
769           IF ( momStepping ) THEN           IF ( momStepping ) THEN
770            CALL CALC_MOM_RHS(            CALL CALC_MOM_RHS(
771       I         bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,       I         bi,bj,iMin,iMax,jMin,jMax,k,km1,kup,kDown,
772       I         xA,yA,uTrans,vTrans,rTrans,rVel,maskC,       I         xA,yA,uTrans,vTrans,rTrans,rVel,maskC,
773       I         phiHyd,KappaRU,KappaRV,       I         phiHyd,KappaRU,KappaRV,
774       U         aTerm,xTerm,cTerm,mTerm,pTerm,       U         aTerm,xTerm,cTerm,mTerm,pTerm,
# Line 669  C--      Calculate accelerations in the Line 784  C--      Calculate accelerations in the
784                 END DO                 END DO
785              END DO              END DO
786  #endif  #endif
787  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
788           ENDIF           ENDIF
789  C--      Calculate active tracer tendencies  C--      Calculate active tracer tendencies
790           IF ( tempStepping ) THEN           IF ( tempStepping ) THEN
791            CALL CALC_GT(            CALL CALC_GT(
792       I         bi,bj,iMin,iMax,jMin,jMax, k,kM1,kUp,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
793       I         xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,       I         xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,
794       I         K13,K23,KappaRT,KapGM,       I         KappaRT,
795       U         aTerm,xTerm,fZon,fMer,fVerT,       U         aTerm,xTerm,fZon,fMer,fVerT,
796       I         myTime, myThid)       I         myTime, myThid)
797           ENDIF           ENDIF
798           IF ( saltStepping ) THEN           IF ( saltStepping ) THEN
799            CALL CALC_GS(            CALL CALC_GS(
800       I         bi,bj,iMin,iMax,jMin,jMax, k,kM1,kUp,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
801       I         xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,       I         xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,
802       I         K13,K23,KappaRS,KapGM,       I         KappaRS,
803       U         aTerm,xTerm,fZon,fMer,fVerS,       U         aTerm,xTerm,fZon,fMer,fVerS,
804       I         myTime, myThid)       I         myTime, myThid)
805           ENDIF           ENDIF
806  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
807  C--      Calculate future values on open boundaries  C--      Calculate future values on open boundaries
808           IF (openBoundaries) THEN           IF (openBoundaries) THEN
809  Caja      CALL CYCLE_OBCS( K, bi, bj, myThid )  Caja      CALL CYCLE_OBCS( k, bi, bj, myThid )
810            CALL SET_OBCS( K, bi, bj, myTime+deltaTclock, myThid )            CALL SET_OBCS( k, bi, bj, myTime+deltaTclock, myThid )
811           ENDIF           ENDIF
812  #endif  #endif
813  C--      Prediction step (step forward all model variables)  C--      Prediction step (step forward all model variables)
814           CALL TIMESTEP(           CALL TIMESTEP(
815       I       bi,bj,iMin,iMax,jMin,jMax,K,       I       bi,bj,iMin,iMax,jMin,jMax,k,
816       I       myIter, myThid)       I       myIter, myThid)
817  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
818  C--      Apply open boundary conditions  C--      Apply open boundary conditions
819           IF (openBoundaries) THEN           IF (openBoundaries) THEN
820  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
821  CADJ STORE gunm1(:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  CADJ STORE gunm1(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
822  CADJ STORE gvnm1(:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  CADJ STORE gvnm1(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
823  CADJ STORE gwnm1(:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  CADJ STORE gwnm1(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
824  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
825              CALL APPLY_OBCS2( bi, bj, K, myThid )  
826                CALL APPLY_OBCS2( bi, bj, k, myThid )
827           END IF           END IF
828  #endif  #endif
829  C--      Freeze water  C--      Freeze water
830           IF (allowFreezing) THEN           IF (allowFreezing) THEN
831  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
832  CADJ STORE gTNm1(:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  CADJ STORE gTNm1(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
833  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
834              CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, K, myThid )              CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )
835           END IF           END IF
836    
837  #ifdef DIVG_IN_DYNAMICS  #ifdef DIVG_IN_DYNAMICS
838  C--      Diagnose barotropic divergence of predicted fields  C--      Diagnose barotropic divergence of predicted fields
839           CALL CALC_DIV_GHAT(           CALL CALC_DIV_GHAT(
840       I       bi,bj,iMin,iMax,jMin,jMax,K,       I       bi,bj,iMin,iMax,jMin,jMax,k,
841       I       xA,yA,       I       xA,yA,
842       I       myThid)       I       myThid)
843  #endif /* DIVG_IN_DYNAMICS */  #endif /* DIVG_IN_DYNAMICS */
# Line 730  C--      Cumulative diagnostic calculati Line 846  C--      Cumulative diagnostic calculati
846  #ifdef INCLUDE_DIAGNOSTICS_INTERFACE_CODE  #ifdef INCLUDE_DIAGNOSTICS_INTERFACE_CODE
847           IF (taveFreq.GT.0.) THEN           IF (taveFreq.GT.0.) THEN
848            CALL DO_TIME_AVERAGES(            CALL DO_TIME_AVERAGES(
849       I                           myTime, myIter, bi, bj, K, kUp, kDown,       I                           myTime, myIter, bi, bj, k, kup, kDown,
850       I                           K13, K23, rVel, KapGM, ConvectCount,       I                           rVel, ConvectCount,
851       I                           myThid )       I                           myThid )
852           ENDIF           ENDIF
853  #endif  #endif
854    
855    
856          ENDDO ! K  C--     k loop
857            ENDDO
 C--     Implicit diffusion  
         IF (implicitDiffusion) THEN  
858    
859  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
860             maximpl = 6             maximpl = 6
861             iikey = ikact*maximpl             iikey = (ikey-1)*maximpl
862  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
863    
864    C--     Implicit diffusion
865            IF (implicitDiffusion) THEN
866    
867           IF (tempStepping) THEN           IF (tempStepping) THEN
868  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
869              idkey = iikey + 1              idkey = iikey + 1
870  #endif  CADJ STORE gTNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
871    #endif /* ALLOW_AUTODIFF_TAMC */
872              CALL IMPLDIFF(              CALL IMPLDIFF(
873       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
874       I         deltaTtracer, KappaRT,recip_HFacC,       I         deltaTtracer, KappaRT,recip_HFacC,
# Line 761  C--     Implicit diffusion Line 879  C--     Implicit diffusion
879           IF (saltStepping) THEN           IF (saltStepping) THEN
880  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
881           idkey = iikey + 2           idkey = iikey + 2
882  #endif  CADJ STORE gSNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
883    #endif /* ALLOW_AUTODIFF_TAMC */
884              CALL IMPLDIFF(              CALL IMPLDIFF(
885       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
886       I         deltaTtracer, KappaRS,recip_HFacC,       I         deltaTtracer, KappaRS,recip_HFacC,
# Line 769  C--     Implicit diffusion Line 888  C--     Implicit diffusion
888       I         myThid )       I         myThid )
889           END IF           END IF
890    
891          ENDIF ! implicitDiffusion  C--     implicitDiffusion
892            ENDIF
893    
894  C--     Implicit viscosity  C--     Implicit viscosity
895          IF (implicitViscosity) THEN          IF (implicitViscosity) THEN
# Line 777  C--     Implicit viscosity Line 897  C--     Implicit viscosity
897           IF (momStepping) THEN           IF (momStepping) THEN
898  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
899           idkey = iikey + 3           idkey = iikey + 3
900  #endif  CADJ STORE gUNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
901    #endif /* ALLOW_AUTODIFF_TAMC */
902            CALL IMPLDIFF(            CALL IMPLDIFF(
903       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
904       I         deltaTmom, KappaRU,recip_HFacW,       I         deltaTmom, KappaRU,recip_HFacW,
# Line 785  C--     Implicit viscosity Line 906  C--     Implicit viscosity
906       I         myThid )       I         myThid )
907  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
908           idkey = iikey + 4           idkey = iikey + 4
909  #endif  CADJ STORE gVNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
910    #endif /* ALLOW_AUTODIFF_TAMC */
911            CALL IMPLDIFF(            CALL IMPLDIFF(
912       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
913       I         deltaTmom, KappaRV,recip_HFacS,       I         deltaTmom, KappaRV,recip_HFacS,
# Line 796  C--     Implicit viscosity Line 918  C--     Implicit viscosity
918    
919  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
920           idkey = iikey + 5           idkey = iikey + 5
921  #endif  CADJ STORE vVelD(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
922    #endif /* ALLOW_AUTODIFF_TAMC */
923            CALL IMPLDIFF(            CALL IMPLDIFF(
924       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
925       I         deltaTmom, KappaRU,recip_HFacW,       I         deltaTmom, KappaRU,recip_HFacW,
# Line 804  C--     Implicit viscosity Line 927  C--     Implicit viscosity
927       I         myThid )       I         myThid )
928  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
929          idkey = iikey + 6          idkey = iikey + 6
930  #endif  CADJ STORE uVelD(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
931    #endif /* ALLOW_AUTODIFF_TAMC */
932            CALL IMPLDIFF(            CALL IMPLDIFF(
933       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
934       I         deltaTmom, KappaRV,recip_HFacS,       I         deltaTmom, KappaRV,recip_HFacS,
# Line 813  C--     Implicit viscosity Line 937  C--     Implicit viscosity
937    
938  #endif  #endif
939    
940           ENDIF ! momStepping  C--      momStepping
941          ENDIF ! implicitViscosity           ENDIF
942    
943    C--     implicitViscosity
944            ENDIF
945    
946         ENDDO         ENDDO
947        ENDDO        ENDDO
948    
 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 )  
   
   
949        RETURN        RETURN
950        END        END
Legend:
Removed from v.1.49  
changed lines
  Added in v.1.55
  ViewVC Help
Powered by ViewVC 1.1.22