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