/[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.54 by heimbach, Mon Nov 13 16:32:57 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 /* 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            ENDIF
442    #endif
443    
444    C--     Start of downward loop
445            DO k=2,Nr
446    
447  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
448           kkey = ikact*(Nr-2+1) + (k-2) + 1           kkey = (ikey-1)*(Nr-2+1) + (k-2) + 1
449  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
450    
451           BOTTOM_LAYER = K .EQ. Nr           BOTTOM_LAYER = k .EQ. Nr
452    
453  #ifdef DO_PIPELINED_CORRECTION_STEP  #ifdef DO_PIPELINED_CORRECTION_STEP
454           IF ( .NOT. BOTTOM_LAYER ) THEN           IF ( .NOT. BOTTOM_LAYER ) THEN
455  C--       Update fields in layer below according to tendency terms  C--       Update fields in layer below according to tendency terms
456            CALL CORRECTION_STEP(            CALL CORRECTION_STEP(
457       I         bi,bj,iMin,iMax,jMin,jMax,K+1,       I         bi,bj,iMin,iMax,jMin,jMax,k+1,
458       I         etaSurfX,etaSurfY,myTime,myThid)       I         etaSurfX,etaSurfY,myTime,myThid)
459  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
460            IF (openBoundaries) THEN            IF (openBoundaries) THEN
461  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
462  CADJ STORE uvel (:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  CADJ STORE uvel (:,:,k,bi,bj) = comlev1_bibj_k
463  CADJ STORE vvel (:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  CADJ &   , key = kkey, byte = isbyte
464  CADJ STORE theta(:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  CADJ STORE vvel (:,:,k,bi,bj) = comlev1_bibj_k
465  CADJ STORE salt(:,:,k,bi,bj)   = comlev1_2d, key = ikey, byte = isbyte  CADJ &   , key = kkey, byte = isbyte
466  #endif  CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k
467               CALL APPLY_OBCS1( bi, bj, K+1, myThid )  CADJ &   , key = kkey, byte = isbyte
468    CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k
469    CADJ &   , key = kkey, byte = isbyte
470    #endif /* ALLOW_AUTODIFF_TAMC */
471                 CALL APPLY_OBCS1( bi, bj, k+1, myThid )
472            END IF            END IF
473  #endif  #endif
474           ENDIF           ENDIF
475  #endif  #endif /* DO_PIPELINED_CORRECTION_STEP */
476    
477  C--      Density of K level (below W(K)) reference to K level  C--      Density of k level (below W(k)) reference to k level
478  #ifdef  INCLUDE_FIND_RHO_CALL  #ifdef  INCLUDE_FIND_RHO_CALL
479  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
480  CADJ STORE theta(:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k
481  CADJ STORE salt (:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  CADJ &   , key = kkey, byte = isbyte
482  #endif  CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k
483    CADJ &   , key = kkey, byte = isbyte
484    #endif /* ALLOW_AUTODIFF_TAMC */
485           CALL FIND_RHO(           CALL FIND_RHO(
486       I      bi, bj, iMin, iMax, jMin, jMax,  K, K, eosType,       I      bi, bj, iMin, iMax, jMin, jMax,  k, k, eosType,
487       O      rhoK,       O      rhoK,
488       I      myThid )       I      myThid )
489    
490    #ifdef ALLOW_AUTODIFF_TAMC
491    cph(   storing not necessary
492    cphCADJ STORE rhoK(:,:)  = comlev1_bibj_k, key = kkey, byte = isbyte
493    cph)
494    #endif /* ALLOW_AUTODIFF_TAMC */
495  #endif  #endif
496           IF (       (.NOT. BOTTOM_LAYER)  
497  #ifdef ALLOW_KPP           IF (.NOT. BOTTOM_LAYER) THEN
498       &       .AND. (.NOT.usingKPPmixing) ! CONVECT not needed with KPP mixing  
 #endif  
      &      ) THEN  
499  C--       Check static stability with layer below and mix as needed.  C--       Check static stability with layer below and mix as needed.
500  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.
501  #ifdef  INCLUDE_FIND_RHO_CALL  #ifdef  INCLUDE_FIND_RHO_CALL
502  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
503  CADJ STORE theta(:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  CADJ STORE theta(:,:,k+1,bi,bj) = comlev1_bibj_k
504  CADJ STORE salt (:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  CADJ &   , key = kkey, byte = isbyte
505  #endif  CADJ STORE salt (:,:,k+1,bi,bj) = comlev1_bibj_k
506    CADJ &   , key = kkey, byte = isbyte
507    #endif /* ALLOW_AUTODIFF_TAMC */
508            CALL FIND_RHO(            CALL FIND_RHO(
509       I       bi, bj, iMin, iMax, jMin, jMax,  K+1, K, eosType,       I       bi, bj, iMin, iMax, jMin, jMax,  k+1, k, eosType,
510       O       rhoKp1,       O       rhoKp1,
511       I       myThid )       I       myThid )
 #endif  
   
512  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
513  CADJ STORE rhok  (:,:)   = comlev1_3d, key = kkey, byte = isbyte  CADJ STORE rhoKp1(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
514  CADJ STORE rhoKm1(:,:)   = comlev1_3d, key = kkey, byte = isbyte  #endif /* ALLOW_AUTODIFF_TAMC */
 CADJ STORE rhoKp1(:,:)   = comlev1_3d, key = kkey, byte = isbyte  
515  #endif  #endif
516    
517  #ifdef  INCLUDE_CONVECT_CALL  #ifdef  INCLUDE_CONVECT_CALL
518            CALL CONVECT(            CALL CONVECT(
519       I        bi,bj,iMin,iMax,jMin,jMax,K+1,rhoK,rhoKp1,       I        bi,bj,iMin,iMax,jMin,jMax,k+1,rhoK,rhoKp1,
520       U        ConvectCount,       U        ConvectCount,
521       I        myTime,myIter,myThid)       I        myTime,myIter,myThid)
522  #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  
523  #endif  #endif
524    
525  C--      Implicit Vertical Diffusion for Convection  C--      Implicit Vertical Diffusion for Convection
526           IF (ivdc_kappa.NE.0.) THEN           IF (ivdc_kappa.NE.0.) THEN
527    #ifdef ALLOW_AUTODIFF_TAMC
528    CADJ STORE rhoKm1(:,:) = comlev1_bibj_k, key = kkey, byte = isbyte
529    #endif /* ALLOW_AUTODIFF_TAMC */
530              CALL CALC_IVDC(              CALL CALC_IVDC(
531       I       bi,bj,iMin,iMax,jMin,jMax,K+1,rhoKm1,rhoKp1,       I       bi,bj,iMin,iMax,jMin,jMax,k+1,rhoKm1,rhoKp1,
532       U       ConvectCount, KappaRT, KappaRS,       U       ConvectCount, KappaRT, KappaRS,
533       I       myTime,myIter,myThid)       I       myTime,myIter,myThid)
 CRG: do we need do store STORE KappaRT, KappaRS ?  
534           END IF           END IF
535    
536  C--       Recompute density after mixing  C--       Recompute density after mixing
537  #ifdef  INCLUDE_FIND_RHO_CALL  #ifdef  INCLUDE_FIND_RHO_CALL
538    #ifdef ALLOW_AUTODIFF_TAMC
539    CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k
540    CADJ &   , key = kkey, byte = isbyte
541    CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k
542    CADJ &   , key = kkey, byte = isbyte
543    #endif /* ALLOW_AUTODIFF_TAMC */
544            CALL FIND_RHO(            CALL FIND_RHO(
545       I       bi, bj, iMin, iMax, jMin, jMax, K, K, eosType,       I       bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,
546       O       rhoK,       O       rhoK,
547       I       myThid )       I       myThid )
548  #endif  #endif
549    
550    C--            IF (.NOT. BOTTOM_LAYER) ends here
551           ENDIF           ENDIF
552    
553  C--      Calculate buoyancy  C--      Calculate buoyancy
554           CALL CALC_BUOYANCY(           CALL CALC_BUOYANCY(
555       I       bi,bj,iMin,iMax,jMin,jMax,K,rhoK,       I       bi,bj,iMin,iMax,jMin,jMax,k,rhoK,
556       O       buoyK,       O       buoyK,
557       I       myThid )       I       myThid )
558    
559  C--      Integrate hydrostatic balance for phiHyd with BC of  C--      Integrate hydrostatic balance for phiHyd with BC of
560  C--      phiHyd(z=0)=0  C--      phiHyd(z=0)=0
561           CALL CALC_PHI_HYD(           CALL CALC_PHI_HYD(
562       I        bi,bj,iMin,iMax,jMin,jMax,K,buoyKm1,buoyK,       I        bi,bj,iMin,iMax,jMin,jMax,k,buoyKm1,buoyK,
563       U        phiHyd,       U        phiHyd,
564       I        myThid )       I        myThid )
565  C--      Calculate iso-neutral slopes for the GM/Redi parameterisation  
566  #ifdef  INCLUDE_FIND_RHO_CALL  #ifdef  INCLUDE_FIND_RHO_CALL
567    C--      Calculate iso-neutral slopes for the GM/Redi parameterisation
568    
569    #ifdef ALLOW_AUTODIFF_TAMC
570    CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k
571    CADJ &   , key = kkey, byte = isbyte
572    CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k
573    CADJ &   , key = kkey, byte = isbyte
574    #endif /* ALLOW_AUTODIFF_TAMC */
575    
576           CALL FIND_RHO(           CALL FIND_RHO(
577       I        bi, bj, iMin, iMax, jMin, jMax, K-1, K, eosType,       I        bi, bj, iMin, iMax, jMin, jMax, k-1, k, eosType,
578       O        rhoTmp,       O        rhoTmp,
579       I        myThid )       I        myThid )
580  #endif  #endif
581    
582  #ifdef  INCLUDE_CALC_ISOSLOPES_CALL  
583  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_GMREDI
584  CADJ STORE rhoTmp(:,:)  = comlev1_3d, key = kkey, byte = isbyte           IF ( useGMRedi ) THEN
585  CADJ STORE rhok  (:,:)  = comlev1_3d, key = kkey, byte = isbyte           CALL GRAD_SIGMA(
586  CADJ STORE rhoKm1(:,:)  = comlev1_3d, key = kkey, byte = isbyte       I             bi, bj, iMin, iMax, jMin, jMax, k,
587  CADJ STORE kapgm (:,:)  = comlev1_3d, key = kkey, byte = isbyte       I             rhoK, rhotmp, rhoK,
588  #endif       O             sigmaX, sigmaY, sigmaR,
589           CALL CALC_ISOSLOPES(       I             myThid )
590       I        bi, bj, iMin, iMax, jMin, jMax, K,           ENDIF
      I        rhoKm1, rhoK, rhotmp,  
      O        K13, K23, K33, KapGM,  
      I        myThid )  
591  #endif  #endif
592    
593           DO J=jMin,jMax           DO J=jMin,jMax
# Line 577  CADJ STORE kapgm (:,:)  = comlev1_3d, ke Line 598  CADJ STORE kapgm (:,:)  = comlev1_3d, ke
598             buoyKm1(I,J) = buoyK(I,J)             buoyKm1(I,J) = buoyK(I,J)
599            ENDDO            ENDDO
600           ENDDO           ENDDO
601          ENDDO  
602  C--     end of k loop  C--     end of k loop
603            ENDDO
604    
605  #ifdef ALLOW_AUTODIFF_TAMC  C     Determines forcing terms based on external fields
606  CADJ STORE theta(:,:,:,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  C     relaxation terms, etc.
607  CADJ STORE salt (:,:,:,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte        CALL EXTERNAL_FORCING_SURF(
608  CADJ STORE uvel (:,:,:,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte       I             bi, bj, iMin, iMax, jMin, jMax,
609  CADJ STORE vvel (:,:,:,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte       I             myThid )
610    
611    #ifdef ALLOW_GMREDI
612            IF (useGMRedi) THEN
613              DO k=1, Nr
614                CALL GMREDI_CALC_TENSOR(
615         I             bi, bj, iMin, iMax, jMin, jMax, k,
616         I             sigmaX, sigmaY, sigmaR,
617         I             myThid )
618              ENDDO
619            ENDIF
620  #endif  #endif
621    
622    #ifdef ALLOW_AUTODIFF_TAMC
623    CADJ STORE KappaRT(:,:,:)     = comlev1_bibj, key = ikey, byte = isbyte
624    CADJ STORE KappaRS(:,:,:)     = comlev1_bibj, key = ikey, byte = isbyte
625    
626    CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
627    CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
628    CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
629    CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
630    
631    C--     dummy initialization to break data flow because
632    C--     calc_div_ghat has a condition for initialization
633            DO J=jMin,jMax
634               DO I=iMin,iMax
635                  cg2d_b(i,j,bi,bj) = 0.0
636               ENDDO
637            ENDDO
638    #endif /* ALLOW_AUTODIFF_TAMC */
639    
640  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
641  C----------------------------------------------  C--   Compute KPP mixing coefficients
642  C--     Compute KPP mixing coefficients          IF (useKPP) THEN
643  C----------------------------------------------  
644          IF (usingKPPmixing) THEN            CALL TIMER_START('KPP_CALC               [DYNAMICS]', myThid)
645  #ifdef ALLOW_AUTODIFF_TAMC            CALL KPP_CALC(
646  CADJ STORE fu  (:,:  ,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte       I                  bi, bj, myTime, myThid )
647  CADJ STORE fv  (:,:  ,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte            CALL TIMER_STOP ('KPP_CALC               [DYNAMICS]', myThid)
648  #endif  
649           CALL TIMER_START('KVMIX (FIND KPP COEFFICIENTS) [DYNAMICS]'  #ifdef ALLOW_AUTODIFF_TAMC
650       I          , myThid)          ELSE
651           CALL KVMIX(            DO j=1-OLy,sNy+OLy
652       I               bi, bj, myTime, myThid )              DO i=1-OLx,sNx+OLx
653           CALL TIMER_STOP ('KVMIX (FIND KPP COEFFICIENTS) [DYNAMICS]'                KPPhbl (i,j,bi,bj) = 1.0
654       I        , myThid)                KPPfrac(i,j,bi,bj) = 0.0
655                  DO k = 1,Nr
656                     KPPghat   (i,j,k,bi,bj) = 0.0
657                     KPPviscAz (i,j,k,bi,bj) = viscAz
658                     KPPdiffKzT(i,j,k,bi,bj) = diffKzT
659                     KPPdiffKzS(i,j,k,bi,bj) = diffKzS
660                  ENDDO
661                ENDDO
662              ENDDO
663    #endif /* ALLOW_AUTODIFF_TAMC */
664          ENDIF          ENDIF
 #endif  
665    
666  C----------------------------------------------  #ifdef ALLOW_AUTODIFF_TAMC
667  C--     start of upward loop  CADJ STORE KPPghat   (:,:,:,bi,bj)
668  C----------------------------------------------  CADJ &   , KPPviscAz (:,:,:,bi,bj)
669          DO K = Nr, 1, -1  CADJ &   , KPPdiffKzT(:,:,:,bi,bj)
670    CADJ &   , KPPdiffKzS(:,:,:,bi,bj)
671           kM1  =max(1,k-1)   ! Points to level above k (=k-1)  CADJ &   , KPPfrac   (:,:  ,bi,bj)
672           kUp  =1+MOD(k+1,2) ! Cycles through 1,2 to point to layer above  CADJ &                 = comlev1_bibj, key = ikey, byte = isbyte
673           kDown=1+MOD(k,2)   ! Cycles through 2,1 to point to current layer  #endif /* ALLOW_AUTODIFF_TAMC */
674    
675    #endif /* ALLOW_KPP */
676    
677    C--     Start of upward loop
678            DO k = Nr, 1, -1
679    
680    C--      km1    Points to level above k (=k-1)
681    C--      kup    Cycles through 1,2 to point to layer above
682    C--      kDown  Cycles through 2,1 to point to current layer
683    
684             km1  =max(1,k-1)
685             kup  =1+MOD(k+1,2)
686             kDown=1+MOD(k,2)
687    
688           iMin = 1-OLx+2           iMin = 1-OLx+2
689           iMax = sNx+OLx-1           iMax = sNx+OLx-1
# Line 620  C--------------------------------------- Line 691  C---------------------------------------
691           jMax = sNy+OLy-1           jMax = sNy+OLy-1
692    
693  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
694           kkey = ikact*(Nr-1+1) + (k-1) + 1           kkey = (ikey-1)*(Nr-1+1) + (k-1) + 1
 #endif  
695    
696  #ifdef ALLOW_AUTODIFF_TAMC  CADJ STORE rvel  (:,:,kDown) = comlev1_bibj_k, key = kkey, byte = isbyte
697  CADJ STORE rvel  (:,:,kDown)  = comlev1_3d, key = kkey, byte = isbyte  CADJ STORE rTrans(:,:)       = comlev1_bibj_k, key = kkey, byte = isbyte
698  CADJ STORE rTrans(:,:)        = comlev1_3d, key = kkey, byte = isbyte  CADJ STORE KappaRT(:,:,:)    = comlev1_bibj_k, key = kkey, byte = isbyte
699  CADJ STORE KappaRT(:,:,:)     = comlev1_3d, key = kkey, byte = isbyte  CADJ STORE KappaRS(:,:,:)    = comlev1_bibj_k, key = kkey, byte = isbyte
700  CADJ STORE KappaRS(:,:,:)     = comlev1_3d, key = kkey, byte = isbyte  #endif /* ALLOW_AUTODIFF_TAMC */
 #endif  
701    
702  C--      Get temporary terms used by tendency routines  C--      Get temporary terms used by tendency routines
703           CALL CALC_COMMON_FACTORS (           CALL CALC_COMMON_FACTORS (
704       I        bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,       I        bi,bj,iMin,iMax,jMin,jMax,k,km1,kup,kDown,
705       O        xA,yA,uTrans,vTrans,rTrans,rVel,maskC,maskUp,       O        xA,yA,uTrans,vTrans,rTrans,rVel,maskC,maskUp,
706       I        myThid)       I        myThid)
707    
708  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
709          IF (openBoundaries) THEN          IF (openBoundaries) THEN
710           CALL APPLY_OBCS3( bi, bj, K, Kup, rTrans, rVel, myThid )           CALL APPLY_OBCS3( bi, bj, k, kup, rTrans, rVel, myThid )
711          ENDIF          ENDIF
712  #endif  #endif
713    
714  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL
715  C--      Calculate the total vertical diffusivity  C--      Calculate the total vertical diffusivity
716           CALL CALC_DIFFUSIVITY(           CALL CALC_DIFFUSIVITY(
717       I        bi,bj,iMin,iMax,jMin,jMax,K,       I        bi,bj,iMin,iMax,jMin,jMax,k,
718       I        maskC,maskUp,KapGM,K33,       I        maskC,maskUp,
719       O        KappaRT,KappaRS,KappaRU,KappaRV,       O        KappaRT,KappaRS,KappaRU,KappaRV,
720       I        myThid)       I        myThid)
721  #endif  #endif
722  C--      Calculate accelerations in the momentum equations  C--      Calculate accelerations in the momentum equations
723           IF ( momStepping ) THEN           IF ( momStepping ) THEN
724            CALL CALC_MOM_RHS(            CALL CALC_MOM_RHS(
725       I         bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,       I         bi,bj,iMin,iMax,jMin,jMax,k,km1,kup,kDown,
726       I         xA,yA,uTrans,vTrans,rTrans,rVel,maskC,       I         xA,yA,uTrans,vTrans,rTrans,rVel,maskC,
727       I         phiHyd,KappaRU,KappaRV,       I         phiHyd,KappaRU,KappaRV,
728       U         aTerm,xTerm,cTerm,mTerm,pTerm,       U         aTerm,xTerm,cTerm,mTerm,pTerm,
# Line 669  C--      Calculate accelerations in the Line 738  C--      Calculate accelerations in the
738                 END DO                 END DO
739              END DO              END DO
740  #endif  #endif
741  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
742           ENDIF           ENDIF
743  C--      Calculate active tracer tendencies  C--      Calculate active tracer tendencies
744           IF ( tempStepping ) THEN           IF ( tempStepping ) THEN
745            CALL CALC_GT(            CALL CALC_GT(
746       I         bi,bj,iMin,iMax,jMin,jMax, k,kM1,kUp,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
747       I         xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,       I         xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,
748       I         K13,K23,KappaRT,KapGM,       I         KappaRT,
749       U         aTerm,xTerm,fZon,fMer,fVerT,       U         aTerm,xTerm,fZon,fMer,fVerT,
750       I         myTime, myThid)       I         myTime, myThid)
751           ENDIF           ENDIF
752           IF ( saltStepping ) THEN           IF ( saltStepping ) THEN
753            CALL CALC_GS(            CALL CALC_GS(
754       I         bi,bj,iMin,iMax,jMin,jMax, k,kM1,kUp,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
755       I         xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,       I         xA,yA,uTrans,vTrans,rTrans,maskUp,maskC,
756       I         K13,K23,KappaRS,KapGM,       I         KappaRS,
757       U         aTerm,xTerm,fZon,fMer,fVerS,       U         aTerm,xTerm,fZon,fMer,fVerS,
758       I         myTime, myThid)       I         myTime, myThid)
759           ENDIF           ENDIF
760  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
761  C--      Calculate future values on open boundaries  C--      Calculate future values on open boundaries
762           IF (openBoundaries) THEN           IF (openBoundaries) THEN
763  Caja      CALL CYCLE_OBCS( K, bi, bj, myThid )  Caja      CALL CYCLE_OBCS( k, bi, bj, myThid )
764            CALL SET_OBCS( K, bi, bj, myTime+deltaTclock, myThid )            CALL SET_OBCS( k, bi, bj, myTime+deltaTclock, myThid )
765           ENDIF           ENDIF
766  #endif  #endif
767  C--      Prediction step (step forward all model variables)  C--      Prediction step (step forward all model variables)
768           CALL TIMESTEP(           CALL TIMESTEP(
769       I       bi,bj,iMin,iMax,jMin,jMax,K,       I       bi,bj,iMin,iMax,jMin,jMax,k,
770       I       myIter, myThid)       I       myIter, myThid)
771  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
772  C--      Apply open boundary conditions  C--      Apply open boundary conditions
773           IF (openBoundaries) THEN           IF (openBoundaries) THEN
774  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
775  CADJ STORE gunm1(:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  CADJ STORE gunm1(:,:,k,bi,bj) = comlev1_bibj_k
776  CADJ STORE gvnm1(:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  CADJ &   , key = kkey, byte = isbyte
777  CADJ STORE gwnm1(:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  CADJ STORE gvnm1(:,:,k,bi,bj) = comlev1_bibj_k
778  #endif  CADJ &   , key = kkey, byte = isbyte
779              CALL APPLY_OBCS2( bi, bj, K, myThid )  CADJ STORE gwnm1(:,:,k,bi,bj) = comlev1_bibj_k
780    CADJ &   , key = kkey, byte = isbyte
781    #endif /* ALLOW_AUTODIFF_TAMC */
782    
783                CALL APPLY_OBCS2( bi, bj, k, myThid )
784           END IF           END IF
785  #endif  #endif
786  C--      Freeze water  C--      Freeze water
787           IF (allowFreezing) THEN           IF (allowFreezing) THEN
788  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
789  CADJ STORE gTNm1(:,:,k,bi,bj)  = comlev1_3d, key = kkey, byte = isbyte  CADJ STORE gTNm1(:,:,k,bi,bj) = comlev1_bibj_k
790  #endif  CADJ &   , key = kkey, byte = isbyte
791              CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, K, myThid )  #endif /* ALLOW_AUTODIFF_TAMC */
792                CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )
793           END IF           END IF
794    
795  #ifdef DIVG_IN_DYNAMICS  #ifdef DIVG_IN_DYNAMICS
796  C--      Diagnose barotropic divergence of predicted fields  C--      Diagnose barotropic divergence of predicted fields
797           CALL CALC_DIV_GHAT(           CALL CALC_DIV_GHAT(
798       I       bi,bj,iMin,iMax,jMin,jMax,K,       I       bi,bj,iMin,iMax,jMin,jMax,k,
799       I       xA,yA,       I       xA,yA,
800       I       myThid)       I       myThid)
801  #endif /* DIVG_IN_DYNAMICS */  #endif /* DIVG_IN_DYNAMICS */
# Line 730  C--      Cumulative diagnostic calculati Line 804  C--      Cumulative diagnostic calculati
804  #ifdef INCLUDE_DIAGNOSTICS_INTERFACE_CODE  #ifdef INCLUDE_DIAGNOSTICS_INTERFACE_CODE
805           IF (taveFreq.GT.0.) THEN           IF (taveFreq.GT.0.) THEN
806            CALL DO_TIME_AVERAGES(            CALL DO_TIME_AVERAGES(
807       I                           myTime, myIter, bi, bj, K, kUp, kDown,       I                           myTime, myIter, bi, bj, k, kup, kDown,
808       I                           K13, K23, rVel, KapGM, ConvectCount,       I                           rVel, ConvectCount,
809       I                           myThid )       I                           myThid )
810           ENDIF           ENDIF
811  #endif  #endif
812    
813    
814          ENDDO ! K  C--     k loop
815            ENDDO
 C--     Implicit diffusion  
         IF (implicitDiffusion) THEN  
816    
817  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
818             maximpl = 6             maximpl = 6
819             iikey = ikact*maximpl             iikey = (ikey-1)*maximpl
820  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
821    
822    C--     Implicit diffusion
823            IF (implicitDiffusion) THEN
824    
825           IF (tempStepping) THEN           IF (tempStepping) THEN
826  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
827              idkey = iikey + 1              idkey = iikey + 1
828  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
829              CALL IMPLDIFF(              CALL IMPLDIFF(
830       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
831       I         deltaTtracer, KappaRT,recip_HFacC,       I         deltaTtracer, KappaRT,recip_HFacC,
# Line 761  C--     Implicit diffusion Line 836  C--     Implicit diffusion
836           IF (saltStepping) THEN           IF (saltStepping) THEN
837  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
838           idkey = iikey + 2           idkey = iikey + 2
839  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
840              CALL IMPLDIFF(              CALL IMPLDIFF(
841       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
842       I         deltaTtracer, KappaRS,recip_HFacC,       I         deltaTtracer, KappaRS,recip_HFacC,
# Line 769  C--     Implicit diffusion Line 844  C--     Implicit diffusion
844       I         myThid )       I         myThid )
845           END IF           END IF
846    
847          ENDIF ! implicitDiffusion  C--     implicitDiffusion
848            ENDIF
849    
850  C--     Implicit viscosity  C--     Implicit viscosity
851          IF (implicitViscosity) THEN          IF (implicitViscosity) THEN
# Line 777  C--     Implicit viscosity Line 853  C--     Implicit viscosity
853           IF (momStepping) THEN           IF (momStepping) THEN
854  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
855           idkey = iikey + 3           idkey = iikey + 3
856  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
857            CALL IMPLDIFF(            CALL IMPLDIFF(
858       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
859       I         deltaTmom, KappaRU,recip_HFacW,       I         deltaTmom, KappaRU,recip_HFacW,
# Line 785  C--     Implicit viscosity Line 861  C--     Implicit viscosity
861       I         myThid )       I         myThid )
862  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
863           idkey = iikey + 4           idkey = iikey + 4
864  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
865            CALL IMPLDIFF(            CALL IMPLDIFF(
866       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
867       I         deltaTmom, KappaRV,recip_HFacS,       I         deltaTmom, KappaRV,recip_HFacS,
# Line 796  C--     Implicit viscosity Line 872  C--     Implicit viscosity
872    
873  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
874           idkey = iikey + 5           idkey = iikey + 5
875  #endif  #endif /* ALLOW_AUTODIFF_TAMC */
876            CALL IMPLDIFF(            CALL IMPLDIFF(
877       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
878       I         deltaTmom, KappaRU,recip_HFacW,       I         deltaTmom, KappaRU,recip_HFacW,
# Line 804  C--     Implicit viscosity Line 880  C--     Implicit viscosity
880       I         myThid )       I         myThid )
881  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
882          idkey = iikey + 6          idkey = iikey + 6
883  #endif  #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         deltaTmom, KappaRV,recip_HFacS,       I         deltaTmom, KappaRV,recip_HFacS,
# Line 813  C--     Implicit viscosity Line 889  C--     Implicit viscosity
889    
890  #endif  #endif
891    
892           ENDIF ! momStepping  C--      momStepping
893          ENDIF ! implicitViscosity           ENDIF
894    
895    C--     implicitViscosity
896            ENDIF
897    
898         ENDDO         ENDDO
899        ENDDO        ENDDO
900    
 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 )  
   
   
901        RETURN        RETURN
902        END        END
Legend:
Removed from v.1.49  
changed lines
  Added in v.1.54
  ViewVC Help
Powered by ViewVC 1.1.22