/[MITgcm]/MITgcm/model/src/thermodynamics.F
ViewVC logotype

Diff of /MITgcm/model/src/thermodynamics.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.68 by heimbach, Fri May 21 21:45:35 2004 UTC revision 1.105 by jmc, Sun Jun 18 23:22:43 2006 UTC
# Line 3  C $Name$ Line 3  C $Name$
3    
4  #include "PACKAGES_CONFIG.h"  #include "PACKAGES_CONFIG.h"
5  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
 #ifdef ALLOW_PTRACERS  
 # include "PTRACERS_OPTIONS.h"  
 #endif  
6    
7  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
8  # ifdef ALLOW_GMREDI  # ifdef ALLOW_GMREDI
# Line 22  C     !INTERFACE: Line 19  C     !INTERFACE:
19        SUBROUTINE THERMODYNAMICS(myTime, myIter, myThid)        SUBROUTINE THERMODYNAMICS(myTime, myIter, myThid)
20  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
21  C     *==========================================================*  C     *==========================================================*
22  C     | SUBROUTINE THERMODYNAMICS                                  C     | SUBROUTINE THERMODYNAMICS
23  C     | o Controlling routine for the prognostic part of the        C     | o Controlling routine for the prognostic part of the
24  C     |   thermo-dynamics.                                          C     |   thermo-dynamics.
25  C     *===========================================================  C     *===========================================================
26  C     | The algorithm...  C     | The algorithm...
27  C     |  C     |
# Line 80  C     == Global variables === Line 77  C     == Global variables ===
77  #include "PARAMS.h"  #include "PARAMS.h"
78  #include "DYNVARS.h"  #include "DYNVARS.h"
79  #include "GRID.h"  #include "GRID.h"
80    #ifdef ALLOW_GENERIC_ADVDIFF
81  #include "GAD.h"  #include "GAD.h"
 #ifdef ALLOW_PASSIVE_TRACER  
 #include "TR1.h"  
82  #endif  #endif
83  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
84    #include "PTRACERS_SIZE.h"
85  #include "PTRACERS.h"  #include "PTRACERS.h"
86  #endif  #endif
87  #ifdef ALLOW_TIMEAVE  #ifdef ALLOW_TIMEAVE
# Line 116  C     myThid - Thread number for this in Line 113  C     myThid - Thread number for this in
113        INTEGER myIter        INTEGER myIter
114        INTEGER myThid        INTEGER myThid
115    
116    #ifdef ALLOW_GENERIC_ADVDIFF
117  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
118  C     == Local variables  C     == Local variables
119  C     xA, yA                 - Per block temporaries holding face areas  C     xA, yA                 - Per block temporaries holding face areas
120  C     uTrans, vTrans, rTrans - Per block temporaries holding flow  C     uFld, vFld, wFld       - Local copy of velocity field (3 components)
121  C                              transport  C     uTrans, vTrans, rTrans - Per block temporaries holding flow transport
122  C                              o uTrans: Zonal transport  C                              o uTrans: Zonal transport
123  C                              o vTrans: Meridional transport  C                              o vTrans: Meridional transport
124  C                              o rTrans: Vertical transport  C                              o rTrans: Vertical transport
# Line 130  C     fVer[STUV]               o fVer: V Line 128  C     fVer[STUV]               o fVer: V
128  C                                      is "pipelined" in the vertical  C                                      is "pipelined" in the vertical
129  C                                      so we need an fVer for each  C                                      so we need an fVer for each
130  C                                      variable.  C                                      variable.
131  C     rhoK, rhoKM1   - Density at current level, and level above  C     kappaRT,       - Total diffusion in vertical at level k, for T and S
132  C     phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean)  C     kappaRS          (background + spatially varying, isopycnal term).
133  C     phiSurfY             or geopotentiel (atmos) in X and Y direction  C     kappaRTr       - Total diffusion in vertical at level k,
134  C     KappaRT,       - Total diffusion in vertical for T and S.  C                      for each passive Tracer
135  C     KappaRS          (background + spatially varying, isopycnal term).  C     kappaRk        - Total diffusion in vertical, all levels, 1 tracer
136  C     useVariableK   = T when vertical diffusion is not constant  C     useVariableK   = T when vertical diffusion is not constant
137  C     iMin, iMax     - Ranges and sub-block indices on which calculations  C     iMin, iMax     - Ranges and sub-block indices on which calculations
138  C     jMin, jMax       are applied.  C     jMin, jMax       are applied.
# Line 144  C     kDown, km1       are switched with Line 142  C     kDown, km1       are switched with
142  C                      index into fVerTerm.  C                      index into fVerTerm.
143        _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
144        _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
145          _RL uFld    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
146          _RL vFld    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
147          _RL wFld    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
148        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
149        _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
150        _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 151  C                      index into fVerTe Line 152  C                      index into fVerTe
152        _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
153        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
154        _RL fVerS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
155  #ifdef ALLOW_PASSIVE_TRACER        _RL kappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
156        _RL fVerTr1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL kappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
 #endif  
157  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
158        _RL fVerP   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2,PTRACERS_num)        _RL fVerP   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2,PTRACERS_num)
159          _RL kappaRTr(1-Olx:sNx+Olx,1-Oly:sNy+Oly,PTRACERS_num)
160  #endif  #endif
161        _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL kappaRk (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
       _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)  
       _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)  
       _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
       _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
       _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
 C     This is currently used by IVDC and Diagnostics  
       _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
       _RL kp1Msk  
       LOGICAL useVariableK  
162        INTEGER iMin, iMax        INTEGER iMin, iMax
163        INTEGER jMin, jMax        INTEGER jMin, jMax
164        INTEGER bi, bj        INTEGER bi, bj
165        INTEGER i, j        INTEGER i, j
166        INTEGER k, km1, kup, kDown        INTEGER k, km1, kup, kDown
167    #ifdef ALLOW_ADAMSBASHFORTH_3
168          INTEGER iterNb, m1, m2
169    #endif
170    #ifdef ALLOW_TIMEAVE
171          LOGICAL useVariableK
172    #endif
173    #ifdef ALLOW_PTRACERS
174        INTEGER iTracer, ip        INTEGER iTracer, ip
175    #endif
176    
177  CEOP  CEOP
178    
# Line 183  CEOP Line 180  CEOP
180           IF ( debugLevel .GE. debLevB )           IF ( debugLevel .GE. debLevB )
181       &    CALL DEBUG_ENTER('THERMODYNAMICS',myThid)       &    CALL DEBUG_ENTER('THERMODYNAMICS',myThid)
182  #endif  #endif
183    
184  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
185  C--   dummy statement to end declaration part  C--   dummy statement to end declaration part
186        ikey = 1        ikey = 1
# Line 201  CHPF$ INDEPENDENT Line 198  CHPF$ INDEPENDENT
198  C--    HPF directive to help TAMC  C--    HPF directive to help TAMC
199  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS
200  CHPF$&                  ,utrans,vtrans,xA,yA  CHPF$&                  ,utrans,vtrans,xA,yA
201  CHPF$&                  ,KappaRT,KappaRS  CHPF$&                  ,kappaRT,kappaRS
202  CHPF$&                  )  CHPF$&                  )
203    # ifdef ALLOW_PTRACERS
204    CHPF$  INDEPENDENT, NEW (fVerP,kappaRTr)
205    # endif
206  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
207    
208         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
# Line 232  C     uninitialised but inert locations. Line 232  C     uninitialised but inert locations.
232            yA(i,j)        = 0. _d 0            yA(i,j)        = 0. _d 0
233            uTrans(i,j)    = 0. _d 0            uTrans(i,j)    = 0. _d 0
234            vTrans(i,j)    = 0. _d 0            vTrans(i,j)    = 0. _d 0
           rhok   (i,j)   = 0. _d 0  
           rhoKM1 (i,j)   = 0. _d 0  
           phiSurfX(i,j)  = 0. _d 0  
           phiSurfY(i,j)  = 0. _d 0  
235            rTrans (i,j)   = 0. _d 0            rTrans (i,j)   = 0. _d 0
236            rTransKp1(i,j) = 0. _d 0            rTransKp1(i,j) = 0. _d 0
237            fVerT  (i,j,1) = 0. _d 0            fVerT  (i,j,1) = 0. _d 0
238            fVerT  (i,j,2) = 0. _d 0            fVerT  (i,j,2) = 0. _d 0
239            fVerS  (i,j,1) = 0. _d 0            fVerS  (i,j,1) = 0. _d 0
240            fVerS  (i,j,2) = 0. _d 0            fVerS  (i,j,2) = 0. _d 0
241  #ifdef ALLOW_PASSIVE_TRACER            kappaRT(i,j)   = 0. _d 0
242            fVerTr1(i,j,1) = 0. _d 0            kappaRS(i,j)   = 0. _d 0
           fVerTr1(i,j,2) = 0. _d 0  
 #endif  
 #ifdef ALLOW_PTRACERS  
           DO ip=1,PTRACERS_num  
            fVerP  (i,j,1,ip) = 0. _d 0  
            fVerP  (i,j,2,ip) = 0. _d 0  
           ENDDO  
 #endif  
243           ENDDO           ENDDO
244          ENDDO          ENDDO
245    
# Line 259  C     uninitialised but inert locations. Line 247  C     uninitialised but inert locations.
247           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
248            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
249  C This is currently also used by IVDC and Diagnostics  C This is currently also used by IVDC and Diagnostics
250             sigmaX(i,j,k) = 0. _d 0             kappaRk(i,j,k)    = 0. _d 0
            sigmaY(i,j,k) = 0. _d 0  
            sigmaR(i,j,k) = 0. _d 0  
            ConvectCount(i,j,k) = 0.  
            KappaRT(i,j,k)    = 0. _d 0  
            KappaRS(i,j,k)    = 0. _d 0  
251  C- tracer tendency needs to be set to zero (moved here from gad_calc_rhs):  C- tracer tendency needs to be set to zero (moved here from gad_calc_rhs):
252             gT(i,j,k,bi,bj)   = 0. _d 0             gT(i,j,k,bi,bj)   = 0. _d 0
253             gS(i,j,k,bi,bj)   = 0. _d 0             gS(i,j,k,bi,bj)   = 0. _d 0
 # ifdef ALLOW_PASSIVE_TRACER  
 ceh3 needs an IF ( use PASSIVE_TRACER) THEN  
            gTr1(i,j,k,bi,bj) = 0. _d 0  
 # endif  
 # ifdef ALLOW_PTRACERS  
 ceh3 this should have an   IF ( usePTRACERS ) THEN  
            DO iTracer=1,PTRACERS_numInUse  
             gPTr(i,j,k,bi,bj,itracer) = 0. _d 0  
            ENDDO  
 # endif  
 #ifdef ALLOW_AUTODIFF_TAMC  
 cph all the following init. are necessary for TAF  
 cph although some of these are re-initialised later.  
 # ifdef ALLOW_GMREDI  
            Kwx(i,j,k,bi,bj)  = 0. _d 0  
            Kwy(i,j,k,bi,bj)  = 0. _d 0  
            Kwz(i,j,k,bi,bj)  = 0. _d 0  
 #  ifdef GM_NON_UNITY_DIAGONAL  
            Kux(i,j,k,bi,bj)  = 0. _d 0  
            Kvy(i,j,k,bi,bj)  = 0. _d 0  
 #  endif  
 #  ifdef GM_EXTRA_DIAGONAL  
            Kuz(i,j,k,bi,bj)  = 0. _d 0  
            Kvz(i,j,k,bi,bj)  = 0. _d 0  
 #  endif  
 #  ifdef GM_BOLUS_ADVEC  
            GM_PsiX(i,j,k,bi,bj)  = 0. _d 0  
            GM_PsiY(i,j,k,bi,bj)  = 0. _d 0  
 #  endif  
 #  ifdef GM_VISBECK_VARIABLE_K  
            VisbeckK(i,j,bi,bj)   = 0. _d 0  
 #  endif  
 # endif /* ALLOW_GMREDI */  
 #endif /* ALLOW_AUTODIFF_TAMC */  
254            ENDDO            ENDDO
255           ENDDO           ENDDO
256          ENDDO          ENDDO
257    
258          iMin = 1-OLx  #ifdef ALLOW_PTRACERS
259          iMax = sNx+OLx          IF ( usePTRACERS ) THEN
260          jMin = 1-OLy           DO ip=1,PTRACERS_num
261          jMax = sNy+OLy             DO j=1-OLy,sNy+OLy
262                DO i=1-OLx,sNx+OLx
263  #ifdef ALLOW_AUTODIFF_TAMC               fVerP  (i,j,1,ip) = 0. _d 0
264  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte               fVerP  (i,j,2,ip) = 0. _d 0
265  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte               kappaRTr(i,j,ip)  = 0. _d 0
266  CADJ STORE totphihyd              ENDDO
267  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte             ENDDO
268  #ifdef ALLOW_KPP           ENDDO
269  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  C-      set tracer tendency to zero:
270  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte           DO iTracer=1,PTRACERS_num
271  #endif            DO k=1,Nr
272  #endif /* ALLOW_AUTODIFF_TAMC */             DO j=1-OLy,sNy+OLy
273                DO i=1-OLx,sNx+OLx
274  #ifdef ALLOW_DEBUG               gPTr(i,j,k,bi,bj,itracer) = 0. _d 0
275          IF ( debugLevel .GE. debLevB )              ENDDO
276       &    CALL DEBUG_MSG('ENTERING UPWARD K LOOP',myThid)             ENDDO
277  #endif            ENDDO
278             ENDDO
279  C--     Start of diagnostic loop          ENDIF
         DO k=Nr,1,-1  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 C? Patrick, is this formula correct now that we change the loop range?  
 C? Do we still need this?  
 cph kkey formula corrected.  
 cph Needed for rhok, rhokm1, in the case useGMREDI.  
          kkey = (itdkey-1)*Nr + k  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 C--       Integrate continuity vertically for vertical velocity  
 c         CALL INTEGRATE_FOR_W(  
 c    I                         bi, bj, k, uVel, vVel,  
 c    O                         wVel,  
 c    I                         myThid )  
   
 #ifdef    ALLOW_OBCS  
 #ifdef    ALLOW_NONHYDROSTATIC  
 C--       Apply OBC to W if in N-H mode  
 c         IF (useOBCS.AND.nonHydrostatic) THEN  
 c           CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )  
 c         ENDIF  
 #endif    /* ALLOW_NONHYDROSTATIC */  
 #endif    /* ALLOW_OBCS */  
   
 C--       Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h  
 C--       MOST of THERMODYNAMICS will be disabled  
 #ifndef SINGLE_LAYER_MODE  
   
 C--       Calculate gradients of potential density for isoneutral  
 C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)  
 c         IF ( k.GT.1 .AND. (useGMRedi.OR.ivdc_kappa.NE.0.) ) THEN  
           IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.) ) THEN  
 #ifdef ALLOW_DEBUG  
             IF ( debugLevel .GE. debLevB )  
      &       CALL DEBUG_CALL('FIND_RHO',myThid)  
 #endif  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  
 CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
             CALL FIND_RHO(  
      I        bi, bj, iMin, iMax, jMin, jMax, k, k,  
      I        theta, salt,  
      O        rhoK,  
      I        myThid )  
   
             IF (k.GT.1) THEN  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  
 CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
              CALL FIND_RHO(  
      I        bi, bj, iMin, iMax, jMin, jMax, k-1, k,  
      I        theta, salt,  
      O        rhoKm1,  
      I        myThid )  
             ENDIF  
 #ifdef ALLOW_DEBUG  
             IF ( debugLevel .GE. debLevB )  
      &       CALL DEBUG_CALL('GRAD_SIGMA',myThid)  
 #endif  
             CALL GRAD_SIGMA(  
      I             bi, bj, iMin, iMax, jMin, jMax, k,  
      I             rhoK, rhoKm1, rhoK,  
      O             sigmaX, sigmaY, sigmaR,  
      I             myThid )  
           ENDIF  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE rhok   (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte  
 CADJ STORE rhokm1 (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
 C--       Implicit Vertical Diffusion for Convection  
 c ==> should use sigmaR !!!  
           IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN  
 #ifdef ALLOW_DEBUG  
             IF ( debugLevel .GE. debLevB )  
      &       CALL DEBUG_CALL('CALC_IVDC',myThid)  
280  #endif  #endif
             CALL CALC_IVDC(  
      I        bi, bj, iMin, iMax, jMin, jMax, k,  
      I        rhoKm1, rhoK,  
      U        ConvectCount, KappaRT, KappaRS,  
      I        myTime, myIter, myThid)  
           ENDIF  
281    
282  #endif /* SINGLE_LAYER_MODE */  #ifdef ALLOW_ADAMSBASHFORTH_3
283    C-      Apply AB on T,S :
284  C--     end of diagnostic k loop (Nr:1)          iterNb = myIter
285          ENDDO          IF (staggerTimeStep) iterNb = myIter - 1
286            m1 = 1 + MOD(iterNb+1,2)
287            m2 = 1 + MOD( iterNb ,2)
288    C       compute T^n+1/2 (stored in gtNm) extrapolating T forward in time
289            IF ( AdamsBashforth_T ) CALL ADAMS_BASHFORTH3(
290         I                                  bi, bj, 0,
291         U                                  theta, gtNm,
292         I                                  tempStartAB, iterNb, myThid )
293    C       compute S^n+1/2 (stored in gsNm) extrapolating S forward in time
294            IF ( AdamsBashforth_S ) CALL ADAMS_BASHFORTH3(
295         I                                  bi, bj, 0,
296         U                                  salt, gsNm,
297         I                                  saltStartAB, iterNb, myThid )
298    #endif /* ALLOW_ADAMSBASHFORTH_3 */
299    
300    c       iMin = 1-OLx
301    c       iMax = sNx+OLx
302    c       jMin = 1-OLy
303    c       jMax = sNy+OLy
304    
305  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
306  cph avoids recomputation of integrate_for_w  cph avoids recomputation of integrate_for_w
307  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
308  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
309    
 #ifdef  ALLOW_OBCS  
 C--     Calculate future values on open boundaries  
         IF (useOBCS) THEN  
 #ifdef ALLOW_DEBUG  
           IF ( debugLevel .GE. debLevB )  
      &     CALL DEBUG_CALL('OBCS_CALC',myThid)  
 #endif  
           CALL OBCS_CALC( bi, bj, myTime+deltaT, myIter+1,  
      I            uVel, vVel, wVel, theta, salt,  
      I            myThid )  
         ENDIF  
 #endif  /* ALLOW_OBCS */  
   
 #ifndef ALLOW_AUTODIFF_TAMC  
         IF ( buoyancyRelation(1:7) .EQ. 'OCEANIC' ) THEN  
 #endif  
 C--     Determines forcing terms based on external fields  
 C       relaxation terms, etc.  
 #ifdef ALLOW_DEBUG  
         IF ( debugLevel .GE. debLevB )  
      &    CALL DEBUG_CALL('EXTERNAL_FORCING_SURF',myThid)  
 #endif  
          CALL EXTERNAL_FORCING_SURF(  
      I             bi, bj, iMin, iMax, jMin, jMax,  
      I             myTime, myIter, myThid )  
 #ifndef ALLOW_AUTODIFF_TAMC  
         ENDIF  
 #endif  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 cph needed for KPP  
 CADJ STORE surfacetendencyU(:,:,bi,bj)  
 CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  
 CADJ STORE surfacetendencyV(:,:,bi,bj)  
 CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  
 CADJ STORE surfacetendencyS(:,:,bi,bj)  
 CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  
 CADJ STORE surfacetendencyT(:,:,bi,bj)  
 CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  
 # ifdef ALLOW_SEAICE  
 CADJ STORE surfacetendencyTice(:,:,bi,bj)  
 CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  
 # endif  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
310  C--     Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h  C--     Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h
311  C--     MOST of THERMODYNAMICS will be disabled  C--     MOST of THERMODYNAMICS will be disabled
312  #ifndef SINGLE_LAYER_MODE  #ifndef SINGLE_LAYER_MODE
313    
 #ifdef  ALLOW_GMREDI  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 cph storing here is needed only for one GMREDI_OPTIONS:  
 cph define GM_BOLUS_ADVEC  
 cph but I've avoided the #ifdef for now, in case more things change  
 CADJ STORE sigmaX(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte  
 CADJ STORE sigmaY(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte  
 CADJ STORE sigmaR(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 C--     Calculate iso-neutral slopes for the GM/Redi parameterisation  
         IF (useGMRedi) THEN  
 #ifdef ALLOW_DEBUG  
           IF ( debugLevel .GE. debLevB )  
      &     CALL DEBUG_CALL('GMREDI_CALC_TENSOR',myThid)  
 #endif  
           CALL GMREDI_CALC_TENSOR(  
      I             bi, bj, iMin, iMax, jMin, jMax,  
      I             sigmaX, sigmaY, sigmaR,  
      I             myThid )  
 #ifdef ALLOW_AUTODIFF_TAMC  
         ELSE  
           CALL GMREDI_CALC_TENSOR_DUMMY(  
      I             bi, bj, iMin, iMax, jMin, jMax,  
      I             sigmaX, sigmaY, sigmaR,  
      I             myThid )  
 #endif /* ALLOW_AUTODIFF_TAMC */  
         ENDIF  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE Kwx(:,:,:,bi,bj)   = comlev1_bibj, key=itdkey, byte=isbyte  
 CADJ STORE Kwy(:,:,:,bi,bj)   = comlev1_bibj, key=itdkey, byte=isbyte  
 CADJ STORE Kwz(:,:,:,bi,bj)   = comlev1_bibj, key=itdkey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 #endif  /* ALLOW_GMREDI */  
   
 #ifdef  ALLOW_KPP  
 C--     Compute KPP mixing coefficients  
         IF (useKPP) THEN  
 #ifdef ALLOW_DEBUG  
           IF ( debugLevel .GE. debLevB )  
      &     CALL DEBUG_CALL('KPP_CALC',myThid)  
 #endif  
           CALL KPP_CALC(  
      I                  bi, bj, myTime, myThid )  
 #ifdef ALLOW_AUTODIFF_TAMC  
         ELSE  
           CALL KPP_CALC_DUMMY(  
      I                  bi, bj, myTime, myThid )  
 #endif /* ALLOW_AUTODIFF_TAMC */  
         ENDIF  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE KPPghat   (:,:,:,bi,bj)  
 CADJ &   , KPPdiffKzT(:,:,:,bi,bj)  
 CADJ &   , KPPdiffKzS(:,:,:,bi,bj)  
 CADJ &   , KPPfrac   (:,:  ,bi,bj)  
 CADJ &                 = comlev1_bibj, key=itdkey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 #endif  /* ALLOW_KPP */  
   
314  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
315  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
316  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
317  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
318  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
319  #ifdef ALLOW_PASSIVE_TRACER  cph-test
320  CADJ STORE tr1  (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  # ifdef ALLOW_DEPTH_CONTROL
321  #endif  CADJ STORE gtnm1(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
322  #ifdef ALLOW_PTRACERS  CADJ STORE gsnm1(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
323  cph-- moved to forward_step to avoid key computation  # endif
 cphCADJ STORE ptracer(:,:,:,bi,bj,itracer) = comlev1_bibj,  
 cphCADJ &                              key=itdkey, byte=isbyte  
 #endif  
324  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
325    
 #ifdef ALLOW_AIM  
 C       AIM - atmospheric intermediate model, physics package code.  
         IF ( useAIM ) THEN  
 #ifdef ALLOW_DEBUG  
           IF ( debugLevel .GE. debLevB )  
      &     CALL DEBUG_CALL('AIM_DO_PHYSICS',myThid)  
 #endif  
          CALL TIMER_START('AIM_DO_PHYSICS   [THERMODYNAMICS]', myThid)  
          CALL AIM_DO_PHYSICS( bi, bj, myTime, myIter, myThid )  
          CALL TIMER_STOP( 'AIM_DO_PHYSICS   [THERMODYNAMICS]', myThid)  
         ENDIF  
 #endif /* ALLOW_AIM */  
   
326  #ifndef DISABLE_MULTIDIM_ADVECTION  #ifndef DISABLE_MULTIDIM_ADVECTION
327  C--     Some advection schemes are better calculated using a multi-dimensional  C--     Some advection schemes are better calculated using a multi-dimensional
328  C       method in the absence of any other terms and, if used, is done here.  C       method in the absence of any other terms and, if used, is done here.
# Line 583  C disable this section of code. Line 340  C disable this section of code.
340       &     CALL DEBUG_CALL('GAD_ADVECTION',myThid)       &     CALL DEBUG_CALL('GAD_ADVECTION',myThid)
341  #endif  #endif
342            CALL GAD_ADVECTION(            CALL GAD_ADVECTION(
343       I             tempImplVertAdv,tempAdvScheme,GAD_TEMPERATURE,       I             tempImplVertAdv, tempAdvScheme, tempVertAdvScheme,
344         I             GAD_TEMPERATURE,
345       I             uVel, vVel, wVel, theta,       I             uVel, vVel, wVel, theta,
346       O             gT,       O             gT,
347       I             bi,bj,myTime,myIter,myThid)       I             bi,bj,myTime,myIter,myThid)
# Line 594  C disable this section of code. Line 352  C disable this section of code.
352       &     CALL DEBUG_CALL('GAD_ADVECTION',myThid)       &     CALL DEBUG_CALL('GAD_ADVECTION',myThid)
353  #endif  #endif
354            CALL GAD_ADVECTION(            CALL GAD_ADVECTION(
355       I             saltImplVertAdv,saltAdvScheme,GAD_SALINITY,       I             saltImplVertAdv, saltAdvScheme, saltVertAdvScheme,
356         I             GAD_SALINITY,
357       I             uVel, vVel, wVel, salt,       I             uVel, vVel, wVel, salt,
358       O             gS,       O             gS,
359       I             bi,bj,myTime,myIter,myThid)       I             bi,bj,myTime,myIter,myThid)
360          ENDIF          ENDIF
361    
362  C Since passive tracers are configurable separately from T,S we  C Since passive tracers are configurable separately from T,S we
363  C call the multi-dimensional method for PTRACERS regardless  C call the multi-dimensional method for PTRACERS regardless
364  C of whether multiDimAdvection is set or not.  C of whether multiDimAdvection is set or not.
# Line 623  C--     Start of thermodynamics loop Line 383  C--     Start of thermodynamics loop
383  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
384  C? Patrick Is this formula correct?  C? Patrick Is this formula correct?
385  cph Yes, but I rewrote it.  cph Yes, but I rewrote it.
386  cph Also, the KappaR? need the index and subscript k!  cph Also, the kappaR? need the index and subscript k!
387           kkey = (itdkey-1)*Nr + k           kkey = (itdkey-1)*Nr + k
388  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
389    
# Line 640  C--       kDown  Cycles through 2,1 to p Line 400  C--       kDown  Cycles through 2,1 to p
400            jMin = 1-OLy            jMin = 1-OLy
401            jMax = sNy+OLy            jMax = sNy+OLy
402    
403            kp1Msk=1.            IF (k.EQ.Nr) THEN
404            IF (k.EQ.Nr) kp1Msk=0.             DO j=1-Oly,sNy+Oly
405                DO i=1-Olx,sNx+Olx
406                 rTransKp1(i,j) = 0. _d 0
407                ENDDO
408               ENDDO
409              ELSE
410             DO j=1-Oly,sNy+Oly             DO j=1-Oly,sNy+Oly
411              DO i=1-Olx,sNx+Olx              DO i=1-Olx,sNx+Olx
412               rTransKp1(i,j) = kp1Msk*rTrans(i,j)               rTransKp1(i,j) = rTrans(i,j)
413              ENDDO              ENDDO
414             ENDDO             ENDDO
415              ENDIF
416  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
417  CADJ STORE rTransKp1(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE rTransKp1(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
418  #endif  #endif
419    
420  C--       Get temporary terms used by tendency routines  C--       Get temporary terms used by tendency routines :
421    C-        Calculate horizontal "volume transport" through tracer cell face
422            CALL CALC_COMMON_FACTORS (            CALL CALC_COMMON_FACTORS (
423       I         bi,bj,iMin,iMax,jMin,jMax,k,       I         uVel, vVel,
424       O         xA,yA,uTrans,vTrans,rTrans,maskUp,       O         uFld, vFld, uTrans, vTrans, xA, yA,
425       I         myThid)       I         k,bi,bj, myThid )
426    
427    C-        Calculate vertical "volume transport" through tracer cell face
428            IF (k.EQ.1) THEN            IF (k.EQ.1) THEN
429  C- Surface interface :  C-         Surface interface :
430             DO j=1-Oly,sNy+Oly             DO j=1-Oly,sNy+Oly
431              DO i=1-Olx,sNx+Olx              DO i=1-Olx,sNx+Olx
432               rTrans(i,j) = 0.               wFld(i,j)   = 0. _d 0
433                 maskUp(i,j) = 0. _d 0
434                 rTrans(i,j) = 0. _d 0
435              ENDDO              ENDDO
436             ENDDO             ENDDO
437            ELSE            ELSE
438  C- Interior interface :  C-         Interior interface :
439             DO j=1-Oly,sNy+Oly             DO j=1-Oly,sNy+Oly
440              DO i=1-Olx,sNx+Olx              DO i=1-Olx,sNx+Olx
441               rTrans(i,j) = rTrans(i,j)*maskC(i,j,k-1,bi,bj)               wFld(i,j)   = wVel(i,j,k,bi,bj)
442                 maskUp(i,j) = maskC(i,j,k-1,bi,bj)*maskC(i,j,k,bi,bj)
443                 rTrans(i,j) = wFld(i,j)*rA(i,j,bi,bj)*maskUp(i,j)
444              ENDDO              ENDDO
445             ENDDO             ENDDO
446            ENDIF            ENDIF
447    
448  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
   
449  C--   Residual transp = Bolus transp + Eulerian transp  C--   Residual transp = Bolus transp + Eulerian transp
450            IF (useGMRedi) THEN            IF (useGMRedi) THEN
451              CALL GMREDI_CALC_UVFLOW(              CALL GMREDI_CALC_UVFLOW(
452       &            uTrans, vTrans, bi, bj, k, myThid)       U                  uFld, vFld, uTrans, vTrans,
453         I                  k, bi, bj, myThid )
454              IF (K.GE.2) CALL GMREDI_CALC_WFLOW(              IF (K.GE.2) CALL GMREDI_CALC_WFLOW(
455       &                    rTrans, bi, bj, k, myThid)       U                  wFld, rTrans,
456         I                  k, bi, bj, myThid )
457            ENDIF            ENDIF
458    # ifdef ALLOW_AUTODIFF_TAMC
 #ifdef ALLOW_AUTODIFF_TAMC  
459  CADJ STORE rTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE rTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
460  #ifdef GM_BOLUS_ADVEC  # ifdef GM_BOLUS_ADVEC
461  CADJ STORE uTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE uTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
462  CADJ STORE vTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE vTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
463  #endif  # endif
464  #endif /* ALLOW_AUTODIFF_TAMC */  # endif /* ALLOW_AUTODIFF_TAMC */
   
465  #endif /* ALLOW_GMREDI */  #endif /* ALLOW_GMREDI */
466    
467  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL
468  C--      Calculate the total vertical diffusivity  C--      Calculate the total vertical diffusivity
469           CALL CALC_DIFFUSIVITY(            IF ( .NOT.implicitDiffusion ) THEN
470       I        bi,bj,iMin,iMax,jMin,jMax,k,              CALL CALC_DIFFUSIVITY(
471       I        maskUp,       I          bi,bj,iMin,iMax,jMin,jMax,k,
472       O        KappaRT,KappaRS,       I          maskUp,
473       I        myThid)       O          kappaRT,kappaRS,
474         I          myThid)
475              ENDIF
476  # ifdef ALLOW_AUTODIFF_TAMC  # ifdef ALLOW_AUTODIFF_TAMC
477  CADJ STORE KappaRT(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE kappaRT(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
478  CADJ STORE KappaRS(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE kappaRS(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
479  # endif /* ALLOW_AUTODIFF_TAMC */  # endif /* ALLOW_AUTODIFF_TAMC */
480  #endif  #endif
481    
# Line 712  CADJ STORE KappaRS(:,:,k)    = comlev1_b Line 485  CADJ STORE KappaRS(:,:,k)    = comlev1_b
485            jMax = sNy+OLy-1            jMax = sNy+OLy-1
486    
487  C--      Calculate active tracer tendencies (gT,gS,...)  C--      Calculate active tracer tendencies (gT,gS,...)
488  C        and step forward storing result in gTnm1, gSnm1, etc.  C        and step forward storing result in gT, gS, etc.
489    C--
490    # ifdef ALLOW_AUTODIFF_TAMC
491    #  if ((defined NONLIN_FRSURF) || (defined ALLOW_DEPTH_CONTROL)) && (defined ALLOW_GMREDI)
492    CADJ STORE kux(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
493    CADJ STORE kvy(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
494    #   ifdef GM_EXTRA_DIAGONAL
495    CADJ STORE kuz(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
496    CADJ STORE kvz(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
497    #   endif
498    #  endif
499    # endif /* ALLOW_AUTODIFF_TAMC */
500           IF ( tempStepping ) THEN           IF ( tempStepping ) THEN
501             CALL CALC_GT(  #ifdef ALLOW_AUTODIFF_TAMC
502    CADJ STORE gTnm1(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
503    # if (defined NONLIN_FRSURF) || (defined ALLOW_DEPTH_CONTROL)
504    CADJ STORE gT(:,:,k,bi,bj)    = comlev1_bibj_k, key=kkey, byte=isbyte
505    # endif
506    #endif /* ALLOW_AUTODIFF_TAMC */
507              CALL CALC_GT(
508       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
509       I         xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,       I         xA, yA, maskUp, uFld, vFld, wFld,
510       I         KappaRT,       I         uTrans, vTrans, rTrans, rTransKp1,
511         I         kappaRT,
512       U         fVerT,       U         fVerT,
513       I         myTime,myIter,myThid)       I         myTime,myIter,myThid)
514    #ifdef ALLOW_ADAMSBASHFORTH_3
515              IF ( AdamsBashforth_T ) THEN
516             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
517       I         bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,       I         bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,
518       I         theta, gT,       I         gtNm(1-Olx,1-Oly,1,1,1,m2),
519         U         gT,
520       I         myIter, myThid)       I         myIter, myThid)
521              ELSE
522    #endif
523               CALL TIMESTEP_TRACER(
524         I         bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,
525         I         theta,
526         U         gT,
527         I         myIter, myThid)
528    #ifdef ALLOW_ADAMSBASHFORTH_3
529              ENDIF
530    #endif
531           ENDIF           ENDIF
532    
533           IF ( saltStepping ) THEN           IF ( saltStepping ) THEN
534             CALL CALC_GS(  #ifdef ALLOW_AUTODIFF_TAMC
535    CADJ STORE gSnm1(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
536    # if (defined NONLIN_FRSURF) || (defined ALLOW_DEPTH_CONTROL)
537    CADJ STORE gS(:,:,k,bi,bj)    = comlev1_bibj_k, key=kkey, byte=isbyte
538    # endif
539    #endif /* ALLOW_AUTODIFF_TAMC */
540              CALL CALC_GS(
541       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
542       I         xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,       I         xA, yA, maskUp, uFld, vFld, wFld,
543       I         KappaRS,       I         uTrans, vTrans, rTrans, rTransKp1,
544         I         kappaRS,
545       U         fVerS,       U         fVerS,
546       I         myTime,myIter,myThid)       I         myTime,myIter,myThid)
547    #ifdef ALLOW_ADAMSBASHFORTH_3
548              IF ( AdamsBashforth_S ) THEN
549             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
550       I         bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,       I         bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,
551       I         salt, gS,       I         gsNm(1-Olx,1-Oly,1,1,1,m2),
552         U         gS,
553       I         myIter, myThid)       I         myIter, myThid)
554           ENDIF            ELSE
555  #ifdef ALLOW_PASSIVE_TRACER  #endif
 ceh3 needs an IF ( usePASSIVE_TRACER ) THEN  
          IF ( tr1Stepping ) THEN  
            CALL CALC_GTR1(  
      I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,  
      I         xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,  
      I         KappaRT,  
      U         fVerTr1,  
      I         myTime,myIter,myThid)  
556             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
557       I         bi,bj,iMin,iMax,jMin,jMax,k,tracerAdvScheme,       I         bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,
558       I         Tr1, gTr1,       I         salt,
559       I         myIter,myThid)       U         gS,
560           ENDIF       I         myIter, myThid)
561    #ifdef ALLOW_ADAMSBASHFORTH_3
562              ENDIF
563  #endif  #endif
564             ENDIF
565    
566  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
567           IF ( usePTRACERS ) THEN           IF ( usePTRACERS ) THEN
568               IF ( .NOT.implicitDiffusion ) THEN
569                 CALL PTRACERS_CALC_DIFF(
570         I            bi,bj,iMin,iMax,jMin,jMax,k,
571         I            maskUp,
572         O            kappaRTr,
573         I            myThid)
574               ENDIF
575    # ifdef ALLOW_AUTODIFF_TAMC
576    CADJ STORE kappaRTr(:,:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
577    # endif /* ALLOW_AUTODIFF_TAMC */
578             CALL PTRACERS_INTEGRATE(             CALL PTRACERS_INTEGRATE(
579       I         bi,bj,k,       I          bi,bj,k,
580       I         xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,       I          xA, yA, maskUp, uFld, vFld, wFld,
581       X         fVerP, KappaRS,       I          uTrans, vTrans, rTrans, rTransKp1,
582       I         myIter,myTime,myThid)       I          kappaRTr,
583         U          fVerP,
584         I          myTime,myIter,myThid)
585           ENDIF           ENDIF
586  #endif /* ALLOW_PTRACERS */  #endif /* ALLOW_PTRACERS */
587    
# Line 785  CADJ &   , key = kkey, byte = isbyte Line 607  CADJ &   , key = kkey, byte = isbyte
607  C--     end of thermodynamic k loop (Nr:1)  C--     end of thermodynamic k loop (Nr:1)
608          ENDDO          ENDDO
609    
610    C       All explicit advection/diffusion/sources should now be
611    C       done. The updated tracer field is in gPtr. Accumalate
612    C       explicit tendency and also reset gPtr to initial tracer
613    C       field for implicit matrix calculation
614    
615    #ifdef ALLOW_MATRIX
616            IF (useMATRIX)
617         &    CALL MATRIX_STORE_TENDENCY_EXP(bi,bj, myTime,myIter,myThid)
618    #endif
619    
620            iMin = 1
621            iMax = sNx
622            jMin = 1
623            jMax = sNy
624    
625  C--     Implicit vertical advection & diffusion  C--     Implicit vertical advection & diffusion
626            IF ( tempStepping .AND. implicitDiffusion ) THEN
627              CALL CALC_3D_DIFFUSIVITY(
628         I         bi,bj,iMin,iMax,jMin,jMax,
629         I         GAD_TEMPERATURE, useGMredi, useKPP,
630         O         kappaRk,
631         I         myThid)
632            ENDIF
633  #ifdef INCLUDE_IMPLVERTADV_CODE  #ifdef INCLUDE_IMPLVERTADV_CODE
634          IF ( tempImplVertAdv ) THEN          IF ( tempImplVertAdv ) THEN
635            CALL GAD_IMPLICIT_R(            CALL GAD_IMPLICIT_R(
636       I         tempImplVertAdv, tempAdvScheme, GAD_TEMPERATURE,       I         tempImplVertAdv, tempAdvScheme, GAD_TEMPERATURE,
637       I         kappaRT, wVel, theta,       I         kappaRk, wVel, theta,
638       U         gT,       U         gT,
639       I         bi, bj, myTime, myIter, myThid )       I         bi, bj, myTime, myIter, myThid )
640          ELSEIF ( tempStepping .AND. implicitDiffusion ) THEN          ELSEIF ( tempStepping .AND. implicitDiffusion ) THEN
# Line 799  C--     Implicit vertical advection & di Line 642  C--     Implicit vertical advection & di
642          IF     ( tempStepping .AND. implicitDiffusion ) THEN          IF     ( tempStepping .AND. implicitDiffusion ) THEN
643  #endif /* INCLUDE_IMPLVERTADV_CODE */  #endif /* INCLUDE_IMPLVERTADV_CODE */
644  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
645  CADJ STORE KappaRT(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte  CADJ STORE kappaRk(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
646  CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte  CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
647  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
648            CALL IMPLDIFF(            CALL IMPLDIFF(
649       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
650       I         deltaTtracer, KappaRT, recip_HFacC,       I         GAD_TEMPERATURE, kappaRk, recip_hFacC,
651       U         gT,       U         gT,
652       I         myThid )       I         myThid )
653          ENDIF          ENDIF
654    
655    #ifdef ALLOW_TIMEAVE
656            useVariableK = useKPP .OR. usePP81 .OR. useMY82 .OR. useGGL90
657         &       .OR. useGMredi .OR. ivdc_kappa.NE.0.
658            IF (taveFreq.GT.0. .AND. useVariableK ) THEN
659             IF (implicitDiffusion) THEN
660               CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, gT, kappaRk,
661         I                        Nr, 3, deltaTclock, bi, bj, myThid)
662    c        ELSE
663    c          CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, theta, kappaRT,
664    c    I                        Nr, 3, deltaTclock, bi, bj, myThid)
665             ENDIF
666            ENDIF
667    #endif /* ALLOW_TIMEAVE */
668    
669            IF ( saltStepping .AND. implicitDiffusion ) THEN
670              CALL CALC_3D_DIFFUSIVITY(
671         I         bi,bj,iMin,iMax,jMin,jMax,
672         I         GAD_SALINITY, useGMredi, useKPP,
673         O         kappaRk,
674         I         myThid)
675            ENDIF
676    
677  #ifdef INCLUDE_IMPLVERTADV_CODE  #ifdef INCLUDE_IMPLVERTADV_CODE
678          IF ( saltImplVertAdv ) THEN          IF ( saltImplVertAdv ) THEN
679            CALL GAD_IMPLICIT_R(            CALL GAD_IMPLICIT_R(
680       I         saltImplVertAdv, saltAdvScheme, GAD_SALINITY,       I         saltImplVertAdv, saltAdvScheme, GAD_SALINITY,
681       I         kappaRS, wVel, salt,       I         kappaRk, wVel, salt,
682       U         gS,       U         gS,
683       I         bi, bj, myTime, myIter, myThid )       I         bi, bj, myTime, myIter, myThid )
684          ELSEIF ( saltStepping .AND. implicitDiffusion ) THEN          ELSEIF ( saltStepping .AND. implicitDiffusion ) THEN
# Line 821  CADJ STORE gT(:,:,:,bi,bj) = comlev1_bib Line 686  CADJ STORE gT(:,:,:,bi,bj) = comlev1_bib
686          IF     ( saltStepping .AND. implicitDiffusion ) THEN          IF     ( saltStepping .AND. implicitDiffusion ) THEN
687  #endif /* INCLUDE_IMPLVERTADV_CODE */  #endif /* INCLUDE_IMPLVERTADV_CODE */
688  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
689  CADJ STORE KappaRS(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte  CADJ STORE kappaRk(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
690  CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte  CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
691  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
692            CALL IMPLDIFF(            CALL IMPLDIFF(
693       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
694       I         deltaTtracer, KappaRS, recip_HFacC,       I         GAD_SALINITY, kappaRk, recip_hFacC,
695       U         gS,       U         gS,
696       I         myThid )       I         myThid )
697          ENDIF          ENDIF
698    
 #ifdef ALLOW_PASSIVE_TRACER  
         IF ( tr1Stepping .AND. implicitDiffusion ) THEN  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE gTr1(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
           CALL IMPLDIFF(  
      I      bi, bj, iMin, iMax, jMin, jMax,  
      I      deltaTtracer, KappaRT, recip_HFacC,  
      U      gTr1,  
      I      myThid )  
         ENDIF  
 #endif  
   
699  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
700  c #ifdef INCLUDE_IMPLVERTADV_CODE          IF     ( usePTRACERS ) THEN
701  c       IF ( usePTRACERS .AND. ptracerImplVertAdv ) THEN  C--     Vertical advection/diffusion (implicit) for passive tracers
702  c       ELSEIF ( usePTRACERS .AND. implicitDiffusion ) THEN             CALL PTRACERS_IMPLICIT(
703  c #else       U                             kappaRk,
704          IF     ( usePTRACERS .AND. implicitDiffusion ) THEN       I                             bi, bj, myTime, myIter, myThid )
 C--     Vertical diffusion (implicit) for passive tracers  
            CALL PTRACERS_IMPLDIFF( bi,bj,KappaRS,myThid )  
705          ENDIF          ENDIF
706  #endif /* ALLOW_PTRACERS */  #endif /* ALLOW_PTRACERS */
707    
# Line 867  C--      Apply open boundary conditions Line 717  C--      Apply open boundary conditions
717          ENDIF          ENDIF
718  #endif   /* ALLOW_OBCS */  #endif   /* ALLOW_OBCS */
719    
 #ifdef ALLOW_TIMEAVE  
 #ifndef HRCUBE  
         IF (taveFreq.GT.0. .AND. ivdc_kappa.NE.0.) THEN  
           CALL TIMEAVE_CUMUL_1T(ConvectCountTave, ConvectCount,  
      I                           Nr, deltaTclock, bi, bj, myThid)  
         ENDIF  
         useVariableK = useKPP .OR. useGMredi .OR. ivdc_kappa.NE.0.  
         IF (taveFreq.GT.0. .AND. useVariableK ) THEN  
          IF (implicitDiffusion) THEN  
           CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, gT, kappaRT,  
      I                        Nr, 3, deltaTclock, bi, bj, myThid)  
          ELSE  
           CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, theta, kappaRT,  
      I                        Nr, 3, deltaTclock, bi, bj, myThid)  
          ENDIF  
         ENDIF  
 #endif /* ndef HRCUBE */  
 #endif /* ALLOW_TIMEAVE */  
   
720  #endif /* SINGLE_LAYER_MODE */  #endif /* SINGLE_LAYER_MODE */
721    
722  C--   end bi,bj loops.  C--   end bi,bj loops.
# Line 899  C--   end bi,bj loops. Line 730  C--   end bi,bj loops.
730         CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (THERMODYNAMICS)',myThid)
731         CALL DEBUG_STATS_RL(Nr,theta,'Theta (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,theta,'Theta (THERMODYNAMICS)',myThid)
732         CALL DEBUG_STATS_RL(Nr,salt,'Salt (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,salt,'Salt (THERMODYNAMICS)',myThid)
733         CALL DEBUG_STATS_RL(Nr,Gt,'Gt (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,gT,'Gt (THERMODYNAMICS)',myThid)
734         CALL DEBUG_STATS_RL(Nr,Gs,'Gs (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,gS,'Gs (THERMODYNAMICS)',myThid)
735         CALL DEBUG_STATS_RL(Nr,GtNm1,'GtNm1 (THERMODYNAMICS)',myThid)  #ifndef ALLOW_ADAMSBASHFORTH_3
736         CALL DEBUG_STATS_RL(Nr,GsNm1,'GsNm1 (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,gtNm1,'GtNm1 (THERMODYNAMICS)',myThid)
737           CALL DEBUG_STATS_RL(Nr,gsNm1,'GsNm1 (THERMODYNAMICS)',myThid)
738    #endif
739  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
740         IF ( usePTRACERS ) THEN         IF ( usePTRACERS ) THEN
741           CALL PTRACERS_DEBUG(myThid)           CALL PTRACERS_DEBUG(myThid)
742         ENDIF         ENDIF
743  #endif /* ALLOW_PTRACERS */  #endif /* ALLOW_PTRACERS */
744        ENDIF        ENDIF
745  #endif  #endif /* ALLOW_DEBUG */
746    
747  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
748           IF ( debugLevel .GE. debLevB )           IF ( debugLevel .GE. debLevB )
749       &    CALL DEBUG_LEAVE('THERMODYNAMICS',myThid)       &    CALL DEBUG_LEAVE('THERMODYNAMICS',myThid)
750  #endif  #endif
751    
752    #endif /* ALLOW_GENERIC_ADVDIFF */
753    
754        RETURN        RETURN
755        END        END
Legend:
Removed from v.1.68  
changed lines
  Added in v.1.105
  ViewVC Help
Powered by ViewVC 1.1.22