/[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.50 by adcroft, Wed Jun 21 19:13:11 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 37  C     == Global variables === Line 31  C     == Global variables ===
31  #include "GRID.h"  #include "GRID.h"
32    
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 90  C     KappaRS          (background + spa Line 87  C     KappaRS          (background + spa
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 129  C                      index into fVerTe Line 126  C                      index into fVerTe
126        _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
127        _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _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 146  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 163  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 198  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 218  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
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
# Line 239  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,  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 270  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 308  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 319  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 355  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       &     ) 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 382  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    #ifdef ALLOW_GMREDI
432          CALL GRAD_SIGMA(          CALL GRAD_SIGMA(
433       I            bi, bj, iMin, iMax, jMin, jMax, K,       I            bi, bj, iMin, iMax, jMin, jMax, k,
434       I            rhoKm1, rhoKm1, rhoKm1,       I            rhoKm1, rhoKm1, rhoKm1,
435       O            sigmaX, sigmaY, sigmaR,       O            sigmaX, sigmaY, sigmaR,
436       I            myThid )       I            myThid )
437    #endif
438    
439  C--     Start of downward loop  C--     Start of downward loop
440          DO K=2,Nr          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       &      ) THEN           IF (.NOT. BOTTOM_LAYER) THEN
491    
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    
576    #ifdef ALLOW_GMREDI
577           CALL GRAD_SIGMA(           CALL GRAD_SIGMA(
578       I             bi, bj, iMin, iMax, jMin, jMax, K,       I             bi, bj, iMin, iMax, jMin, jMax, k,
579       I             rhoK, rhotmp, rhoK,       I             rhoK, rhotmp, rhoK,
580       O             sigmaX, sigmaY, sigmaR,       O             sigmaX, sigmaY, sigmaR,
581       I             myThid )       I             myThid )
582    #endif
583    
584           DO J=jMin,jMax           DO J=jMin,jMax
585            DO I=iMin,iMax            DO I=iMin,iMax
# Line 554  C--      Calculate iso-neutral slopes fo Line 589  C--      Calculate iso-neutral slopes fo
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_GMREDI  #ifdef ALLOW_GMREDI
597  #ifdef ALLOW_AUTODIFF_TAMC          IF (useGMRedi) THEN
598  CADJ STORE rhoTmp(:,:)  = comlev1_3d, key = kkey, byte = isbyte            DO k=1, Nr
599  CADJ STORE rhok  (:,:)  = comlev1_3d, key = kkey, byte = isbyte              CALL GMREDI_CALC_TENSOR(
600  CADJ STORE rhoKm1(:,:)  = comlev1_3d, key = kkey, byte = isbyte       I             bi, bj, iMin, iMax, jMin, jMax, k,
 #endif  
         DO K=1, Nr  
          IF (use_GMRedi) CALL GMREDI_CALC_TENSOR(  
      I             bi, bj, iMin, iMax, jMin, jMax, K,  
601       I             sigmaX, sigmaY, sigmaR,       I             sigmaX, sigmaY, sigmaR,
602       I             myThid )       I             myThid )
603          ENDDO            ENDDO
604            ENDIF
605  #endif  #endif
606    
607  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
608  CADJ STORE theta(:,:,:,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  CADJ STORE KappaRT(:,:,:)     = comlev1_bibj, key = ikey, byte = isbyte
609  CADJ STORE salt (:,:,:,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  CADJ STORE KappaRS(:,:,:)     = comlev1_bibj, key = ikey, byte = isbyte
610  CADJ STORE uvel (:,:,:,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  
611  CADJ STORE vvel (:,:,:,bi,bj)  = comlev1_2d, key = ikey, byte = isbyte  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte
612  #endif  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--     Compute KPP mixing coefficients  C--   Compute KPP mixing coefficients
628          CALL TIMER_START('KPP_CALC               [DYNAMICS]', myThid)          IF (useKPP) THEN
629          CALL KPP_CALC(  
630       I               bi, bj, myTime, myThid )            CALL TIMER_START('KPP_CALC               [DYNAMICS]', myThid)
631          CALL TIMER_STOP ('KPP_CALC               [DYNAMICS]', myThid)            CALL KPP_CALC(
632  #endif       I                 bi, bj, myTime, myThid )
633              CALL TIMER_STOP ('KPP_CALC               [DYNAMICS]', myThid)
634    
635    #ifdef ALLOW_AUTODIFF_TAMC
636            ELSE
637              DO j=1-OLy,sNy+OLy
638                DO i=1-OLx,sNx+OLx
639                  KPPhbl (i,j,bi,bj) = 1.0
640                  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
651    
652    #ifdef ALLOW_AUTODIFF_TAMC
653    CADJ STORE KPPghat   (:,:,:,bi,bj)
654    CADJ &   , KPPviscAz (:,:,:,bi,bj)
655    CADJ &   , KPPdiffKzT(:,:,:,bi,bj)
656    CADJ &   , KPPdiffKzS(:,:,:,bi,bj)
657    CADJ &   , KPPfrac   (:,:  ,bi,bj)
658    CADJ &                 = comlev1_bibj, key = ikey, byte = isbyte
659    #endif /* ALLOW_AUTODIFF_TAMC */
660    
661    #endif /* ALLOW_KPP */
662    
663  C--     Start of upward loop  C--     Start of upward loop
664          DO K = Nr, 1, -1          DO k = Nr, 1, -1
665    
666           kM1  =max(1,k-1)   ! Points to level above k (=k-1)  C--      km1    Points to level above k (=k-1)
667           kUp  =1+MOD(k+1,2) ! Cycles through 1,2 to point to layer above  C--      kup    Cycles through 1,2 to point to layer above
668           kDown=1+MOD(k,2)   ! Cycles through 2,1 to point to current layer  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 599  C--     Start of upward loop Line 677  C--     Start of upward loop
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,       I        maskC,maskUp,
705       O        KappaRT,KappaRS,KappaRU,KappaRV,       O        KappaRT,KappaRS,KappaRU,KappaRV,
706       I        myThid)       I        myThid)
# Line 632  C--      Calculate the total vertical di Line 708  C--      Calculate the total vertical di
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 648  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         KappaRT,       I         KappaRT,
735       U         aTerm,xTerm,fZon,fMer,fVerT,       U         aTerm,xTerm,fZon,fMer,fVerT,
# Line 661  C--      Calculate active tracer tendenc Line 737  C--      Calculate active tracer tendenc
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         KappaRS,       I         KappaRS,
743       U         aTerm,xTerm,fZon,fMer,fVerS,       U         aTerm,xTerm,fZon,fMer,fVerS,
# Line 670  C--      Calculate active tracer tendenc Line 746  C--      Calculate active tracer tendenc
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 709  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                           rVel, 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 740  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 748  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 756  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 764  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 775  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 783  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 792  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.)  
 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.50  
changed lines
  Added in v.1.53
  ViewVC Help
Powered by ViewVC 1.1.22