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