/[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.65 by dimitri, Sun Jan 25 00:31:52 2004 UTC revision 1.99 by heimbach, Mon Mar 6 18:25:49 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 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"
82  #ifdef ALLOW_PASSIVE_TRACER  #endif
83  #include "TR1.h"  #ifdef ALLOW_OFFLINE
84    #include "OFFLINE.h"
85  #endif  #endif
86  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
87    #include "PTRACERS_SIZE.h"
88  #include "PTRACERS.h"  #include "PTRACERS.h"
89  #endif  #endif
90  #ifdef ALLOW_TIMEAVE  #ifdef ALLOW_TIMEAVE
# Line 102  C     == Global variables === Line 102  C     == Global variables ===
102  # ifdef ALLOW_GMREDI  # ifdef ALLOW_GMREDI
103  #  include "GMREDI.h"  #  include "GMREDI.h"
104  # endif  # endif
105    # ifdef ALLOW_EBM
106    #  include "EBM.h"
107    # endif
108  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
109    
110  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
# Line 113  C     myThid - Thread number for this in Line 116  C     myThid - Thread number for this in
116        INTEGER myIter        INTEGER myIter
117        INTEGER myThid        INTEGER myThid
118    
119    #ifdef ALLOW_GENERIC_ADVDIFF
120  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
121  C     == Local variables  C     == Local variables
122  C     xA, yA                 - Per block temporaries holding face areas  C     xA, yA                 - Per block temporaries holding face areas
# Line 127  C     fVer[STUV]               o fVer: V Line 131  C     fVer[STUV]               o fVer: V
131  C                                      is "pipelined" in the vertical  C                                      is "pipelined" in the vertical
132  C                                      so we need an fVer for each  C                                      so we need an fVer for each
133  C                                      variable.  C                                      variable.
134  C     rhoK, rhoKM1   - Density at current level, and level above  C     kappaRT,       - Total diffusion in vertical at level k, for T and S
135  C     phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean)  C     kappaRS          (background + spatially varying, isopycnal term).
136  C     phiSurfY             or geopotentiel (atmos) in X and Y direction  C     kappaRTr       - Total diffusion in vertical at level k,
137  C     KappaRT,       - Total diffusion in vertical for T and S.  C                      for each passive Tracer
138  C     KappaRS          (background + spatially varying, isopycnal term).  C     kappaRk        - Total diffusion in vertical, all levels, 1 tracer  
139  C     useVariableK   = T when vertical diffusion is not constant  C     useVariableK   = T when vertical diffusion is not constant
140  C     iMin, iMax     - Ranges and sub-block indices on which calculations  C     iMin, iMax     - Ranges and sub-block indices on which calculations
141  C     jMin, jMax       are applied.  C     jMin, jMax       are applied.
# Line 148  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 198  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  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
204    
# Line 229  C     uninitialised but inert locations. Line 229  C     uninitialised but inert locations.
229            yA(i,j)        = 0. _d 0            yA(i,j)        = 0. _d 0
230            uTrans(i,j)    = 0. _d 0            uTrans(i,j)    = 0. _d 0
231            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  
232            rTrans (i,j)   = 0. _d 0            rTrans (i,j)   = 0. _d 0
233            rTransKp1(i,j) = 0. _d 0            rTransKp1(i,j) = 0. _d 0
234            fVerT  (i,j,1) = 0. _d 0            fVerT  (i,j,1) = 0. _d 0
235            fVerT  (i,j,2) = 0. _d 0            fVerT  (i,j,2) = 0. _d 0
236            fVerS  (i,j,1) = 0. _d 0            fVerS  (i,j,1) = 0. _d 0
237            fVerS  (i,j,2) = 0. _d 0            fVerS  (i,j,2) = 0. _d 0
238  #ifdef ALLOW_PASSIVE_TRACER            kappaRT(i,j)   = 0. _d 0
239            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  
240           ENDDO           ENDDO
241          ENDDO          ENDDO
242    
# Line 256  C     uninitialised but inert locations. Line 244  C     uninitialised but inert locations.
244           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
245            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
246  C This is currently also used by IVDC and Diagnostics  C This is currently also used by IVDC and Diagnostics
247             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  
248  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):
249             gT(i,j,k,bi,bj)   = 0. _d 0             gT(i,j,k,bi,bj)   = 0. _d 0
250             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 */  
251            ENDDO            ENDDO
252           ENDDO           ENDDO
253          ENDDO          ENDDO
254    
255          iMin = 1-OLx  #ifdef ALLOW_PTRACERS
256          iMax = sNx+OLx          IF ( usePTRACERS ) THEN
257          jMin = 1-OLy           DO ip=1,PTRACERS_num
258          jMax = sNy+OLy             DO j=1-OLy,sNy+OLy
259                DO i=1-OLx,sNx+OLx
260  #ifdef ALLOW_AUTODIFF_TAMC               fVerP  (i,j,1,ip) = 0. _d 0
261  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte               fVerP  (i,j,2,ip) = 0. _d 0
262  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte               kappaRTr(i,j,ip)  = 0. _d 0
263  CADJ STORE totphihyd              ENDDO
264  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte             ENDDO
265  #ifdef ALLOW_KPP           ENDDO
266  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  C-      set tracer tendency to zero:
267  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte           DO iTracer=1,PTRACERS_numInUse
268  #endif            DO k=1,Nr
269  #endif /* ALLOW_AUTODIFF_TAMC */             DO j=1-OLy,sNy+OLy
270                DO i=1-OLx,sNx+OLx
271  #ifdef ALLOW_DEBUG               gPTr(i,j,k,bi,bj,itracer) = 0. _d 0
272          IF ( debugLevel .GE. debLevB )              ENDDO
273       &    CALL DEBUG_MSG('ENTERING UPWARD K LOOP',myThid)             ENDDO
274  #endif            ENDDO
275             ENDDO
276  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)  
277  #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  
   
 #endif /* SINGLE_LAYER_MODE */  
278    
279  C--     end of diagnostic k loop (Nr:1)  #ifdef ALLOW_ADAMSBASHFORTH_3
280          ENDDO  C-      Apply AB on T,S :
281            iterNb = myIter
282            IF (staggerTimeStep) iterNb = myIter - 1
283            m1 = 1 + MOD(iterNb+1,2)
284            m2 = 1 + MOD( iterNb ,2)
285    C       compute T^n+1/2 (stored in gtNm) extrapolating T forward in time
286            IF ( AdamsBashforth_T ) CALL ADAMS_BASHFORTH3(
287         I                                  bi, bj, 0,
288         U                                  theta, gtNm,
289         I                                  tempStartAB, iterNb, myThid )
290    C       compute S^n+1/2 (stored in gsNm) extrapolating S forward in time
291            IF ( AdamsBashforth_S ) CALL ADAMS_BASHFORTH3(
292         I                                  bi, bj, 0,
293         U                                  salt, gsNm,
294         I                                  saltStartAB, iterNb, myThid )
295    #endif /* ALLOW_ADAMSBASHFORTH_3 */
296    
297    c       iMin = 1-OLx
298    c       iMax = sNx+OLx
299    c       jMin = 1-OLy
300    c       jMax = sNy+OLy
301    
302  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
303  cph avoids recomputation of integrate_for_w  cph avoids recomputation of integrate_for_w
304  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
305  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
306    
 #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 */  
   
         IF ( buoyancyRelation(1:7) .EQ. 'OCEANIC' ) THEN  
 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 )  
   
 #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 */  
         ENDIF  
   
307  C--     Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h  C--     Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h
308  C--     MOST of THERMODYNAMICS will be disabled  C--     MOST of THERMODYNAMICS will be disabled
309  #ifndef SINGLE_LAYER_MODE  #ifndef SINGLE_LAYER_MODE
310    
 #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 */  
   
311  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
312  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
313  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
314  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
315  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
316  #ifdef ALLOW_PASSIVE_TRACER  cph-test
317  CADJ STORE tr1  (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  cphCADJ STORE gtnm1(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
318  #endif  cphCADJ STORE gsnm1(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
319  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
320  cph-- moved to forward_step to avoid key computation  cph-- moved to forward_step to avoid key computation
321  cphCADJ STORE ptracer(:,:,:,bi,bj,itracer) = comlev1_bibj,  cphCADJ STORE ptracer(:,:,:,bi,bj,itracer) = comlev1_bibj,
# Line 546  cphCADJ &                              k Line 323  cphCADJ &                              k
323  #endif  #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 570  C to be able to exclude this scheme to a Line 334  C to be able to exclude this scheme to a
334  C recomputation. It *is* differentiable, if you need it.  C recomputation. It *is* differentiable, if you need it.
335  C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to  C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to
336  C disable this section of code.  C disable this section of code.
337    #ifndef ALLOW_OFFLINE
338          IF (tempMultiDimAdvec) THEN          IF (tempMultiDimAdvec) THEN
339  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
340            IF ( debugLevel .GE. debLevB )            IF ( debugLevel .GE. debLevB )
341       &     CALL DEBUG_CALL('GAD_ADVECTION',myThid)       &     CALL DEBUG_CALL('GAD_ADVECTION',myThid)
342  #endif  #endif
343            CALL GAD_ADVECTION(            CALL GAD_ADVECTION(
344       I             tempImplVertAdv,tempAdvScheme,GAD_TEMPERATURE,       I             tempImplVertAdv, tempAdvScheme, tempVertAdvScheme,
345         I             GAD_TEMPERATURE,
346       I             uVel, vVel, wVel, theta,       I             uVel, vVel, wVel, theta,
347       O             gT,       O             gT,
348       I             bi,bj,myTime,myIter,myThid)       I             bi,bj,myTime,myIter,myThid)
349          ENDIF          ENDIF
350    #endif
351    #ifndef ALLOW_OFFLINE
352          IF (saltMultiDimAdvec) THEN          IF (saltMultiDimAdvec) THEN
353  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
354            IF ( debugLevel .GE. debLevB )            IF ( debugLevel .GE. debLevB )
355       &     CALL DEBUG_CALL('GAD_ADVECTION',myThid)       &     CALL DEBUG_CALL('GAD_ADVECTION',myThid)
356  #endif  #endif
357            CALL GAD_ADVECTION(            CALL GAD_ADVECTION(
358       I             saltImplVertAdv,saltAdvScheme,GAD_SALINITY,       I             saltImplVertAdv, saltAdvScheme, saltVertAdvScheme,
359         I             GAD_SALINITY,
360       I             uVel, vVel, wVel, salt,       I             uVel, vVel, wVel, salt,
361       O             gS,       O             gS,
362       I             bi,bj,myTime,myIter,myThid)       I             bi,bj,myTime,myIter,myThid)
363          ENDIF          ENDIF
364    #endif
365  C Since passive tracers are configurable separately from T,S we  C Since passive tracers are configurable separately from T,S we
366  C call the multi-dimensional method for PTRACERS regardless  C call the multi-dimensional method for PTRACERS regardless
367  C of whether multiDimAdvection is set or not.  C of whether multiDimAdvection is set or not.
# Line 616  C--     Start of thermodynamics loop Line 386  C--     Start of thermodynamics loop
386  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
387  C? Patrick Is this formula correct?  C? Patrick Is this formula correct?
388  cph Yes, but I rewrote it.  cph Yes, but I rewrote it.
389  cph Also, the KappaR? need the index and subscript k!  cph Also, the kappaR? need the index and subscript k!
390           kkey = (itdkey-1)*Nr + k           kkey = (itdkey-1)*Nr + k
391  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
392    
# Line 633  C--       kDown  Cycles through 2,1 to p Line 403  C--       kDown  Cycles through 2,1 to p
403            jMin = 1-OLy            jMin = 1-OLy
404            jMax = sNy+OLy            jMax = sNy+OLy
405    
406            kp1Msk=1.            IF (k.EQ.Nr) THEN
           IF (k.EQ.Nr) kp1Msk=0.  
407             DO j=1-Oly,sNy+Oly             DO j=1-Oly,sNy+Oly
408              DO i=1-Olx,sNx+Olx              DO i=1-Olx,sNx+Olx
409               rTransKp1(i,j) = kp1Msk*rTrans(i,j)               rTransKp1(i,j) = 0. _d 0
410              ENDDO              ENDDO
411             ENDDO             ENDDO
412              ELSE
413               DO j=1-Oly,sNy+Oly
414                DO i=1-Olx,sNx+Olx
415                 rTransKp1(i,j) = rTrans(i,j)
416                ENDDO
417               ENDDO
418              ENDIF
419    #ifdef ALLOW_AUTODIFF_TAMC
420    CADJ STORE rTransKp1(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
421    #endif
422    
423  C--       Get temporary terms used by tendency routines  C--       Get temporary terms used by tendency routines
424            CALL CALC_COMMON_FACTORS (            CALL CALC_COMMON_FACTORS (
# Line 673  C--   Residual transp = Bolus transp + E Line 452  C--   Residual transp = Bolus transp + E
452       &                    rTrans, bi, bj, k, myThid)       &                    rTrans, bi, bj, k, myThid)
453            ENDIF            ENDIF
454    
455  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
456    CADJ STORE rTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
457  #ifdef GM_BOLUS_ADVEC  #ifdef GM_BOLUS_ADVEC
458  CADJ STORE uTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE uTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
459  CADJ STORE vTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE vTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
 CADJ STORE rTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte  
460  #endif  #endif
461  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
462    
# Line 685  CADJ STORE rTrans(:,:)    = comlev1_bibj Line 464  CADJ STORE rTrans(:,:)    = comlev1_bibj
464    
465  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL
466  C--      Calculate the total vertical diffusivity  C--      Calculate the total vertical diffusivity
467           CALL CALC_DIFFUSIVITY(            IF ( .NOT.implicitDiffusion ) THEN
468       I        bi,bj,iMin,iMax,jMin,jMax,k,              CALL CALC_DIFFUSIVITY(
469       I        maskUp,       I          bi,bj,iMin,iMax,jMin,jMax,k,
470       O        KappaRT,KappaRS,       I          maskUp,
471       I        myThid)       O          kappaRT,kappaRS,
472         I          myThid)
473              ENDIF
474  # ifdef ALLOW_AUTODIFF_TAMC  # ifdef ALLOW_AUTODIFF_TAMC
475  CADJ STORE KappaRT(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE kappaRT(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
476  CADJ STORE KappaRS(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE kappaRS(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
477  # endif /* ALLOW_AUTODIFF_TAMC */  # endif /* ALLOW_AUTODIFF_TAMC */
478  #endif  #endif
479    
# Line 702  CADJ STORE KappaRS(:,:,k)    = comlev1_b Line 483  CADJ STORE KappaRS(:,:,k)    = comlev1_b
483            jMax = sNy+OLy-1            jMax = sNy+OLy-1
484    
485  C--      Calculate active tracer tendencies (gT,gS,...)  C--      Calculate active tracer tendencies (gT,gS,...)
486  C        and step forward storing result in gTnm1, gSnm1, etc.  C        and step forward storing result in gT, gS, etc.
487    #ifndef ALLOW_OFFLINE
488           IF ( tempStepping ) THEN           IF ( tempStepping ) THEN
489    #ifdef ALLOW_AUTODIFF_TAMC
490    CADJ STORE gTnm1(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
491    # ifdef NONLIN_FRSURF
492    CADJ STORE gT(:,:,k,bi,bj)    = comlev1_bibj_k, key=kkey, byte=isbyte
493    # endif
494    #endif /* ALLOW_AUTODIFF_TAMC */
495             CALL CALC_GT(             CALL CALC_GT(
496       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
497       I         xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,       I         xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,
498       I         KappaRT,       I         kappaRT,
499       U         fVerT,       U         fVerT,
500       I         myTime,myIter,myThid)       I         myTime,myIter,myThid)
501    #ifdef ALLOW_ADAMSBASHFORTH_3
502              IF ( AdamsBashforth_T ) THEN
503             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
504       I         bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,       I         bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,
505       I         theta, gT,       I         gtNm(1-Olx,1-Oly,1,1,1,m2),
506         U         gT,
507       I         myIter, myThid)       I         myIter, myThid)
508              ELSE
509    #endif
510               CALL TIMESTEP_TRACER(
511         I         bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,
512         I         theta,
513         U         gT,
514         I         myIter, myThid)
515    #ifdef ALLOW_ADAMSBASHFORTH_3
516              ENDIF
517    #endif
518           ENDIF           ENDIF
519    #endif
520    
521    #ifndef ALLOW_OFFLINE
522           IF ( saltStepping ) THEN           IF ( saltStepping ) THEN
523    #ifdef ALLOW_AUTODIFF_TAMC
524    CADJ STORE gSnm1(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
525    # ifdef NONLIN_FRSURF
526    CADJ STORE gS(:,:,k,bi,bj)    = comlev1_bibj_k, key=kkey, byte=isbyte
527    # endif
528    #endif /* ALLOW_AUTODIFF_TAMC */
529             CALL CALC_GS(             CALL CALC_GS(
530       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
531       I         xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,       I         xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,
532       I         KappaRS,       I         kappaRS,
533       U         fVerS,       U         fVerS,
534       I         myTime,myIter,myThid)       I         myTime,myIter,myThid)
535    #ifdef ALLOW_ADAMSBASHFORTH_3
536              IF ( AdamsBashforth_S ) THEN
537             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
538       I         bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,       I         bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,
539       I         salt, gS,       I         gsNm(1-Olx,1-Oly,1,1,1,m2),
540         U         gS,
541       I         myIter, myThid)       I         myIter, myThid)
542           ENDIF            ELSE
543  #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)  
544             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
545       I         bi,bj,iMin,iMax,jMin,jMax,k,tracerAdvScheme,       I         bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,
546       I         Tr1, gTr1,       I         salt,
547       I         myIter,myThid)       U         gS,
548         I         myIter, myThid)
549    #ifdef ALLOW_ADAMSBASHFORTH_3
550              ENDIF
551    #endif
552           ENDIF           ENDIF
553  #endif  #endif
554  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
555           IF ( usePTRACERS ) THEN           IF ( usePTRACERS ) THEN
556               IF ( .NOT.implicitDiffusion ) THEN
557                 CALL PTRACERS_CALC_DIFF(
558         I            bi,bj,iMin,iMax,jMin,jMax,k,
559         I            maskUp,
560         O            kappaRTr,
561         I            myThid)
562               ENDIF
563    # ifdef ALLOW_AUTODIFF_TAMC
564    CADJ STORE kappaRTr(:,:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
565    # endif /* ALLOW_AUTODIFF_TAMC */
566             CALL PTRACERS_INTEGRATE(             CALL PTRACERS_INTEGRATE(
567       I         bi,bj,k,       I         bi,bj,k,
568       I         xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,       I         xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,
569       X         fVerP, KappaRS,       U         fVerP,
570         I         kappaRTr,
571       I         myIter,myTime,myThid)       I         myIter,myTime,myThid)
572           ENDIF           ENDIF
573  #endif /* ALLOW_PTRACERS */  #endif /* ALLOW_PTRACERS */
# Line 763  C--      Apply open boundary conditions Line 582  C--      Apply open boundary conditions
582  C--      Freeze water  C--      Freeze water
583  C  this bit of code is left here for backward compatibility.  C  this bit of code is left here for backward compatibility.
584  C  freezing at surface level has been moved to FORWARD_STEP  C  freezing at surface level has been moved to FORWARD_STEP
585    #ifndef ALLOW_OFFLINE
586           IF ( useOldFreezing .AND. .NOT. useSEAICE           IF ( useOldFreezing .AND. .NOT. useSEAICE
587       &       .AND. .NOT.(useThSIce.AND.k.EQ.1) ) THEN       &       .AND. .NOT.(useThSIce.AND.k.EQ.1) ) THEN
588  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 771  CADJ &   , key = kkey, byte = isbyte Line 591  CADJ &   , key = kkey, byte = isbyte
591  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
592              CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )              CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )
593           ENDIF           ENDIF
594    #endif
595    
596  C--     end of thermodynamic k loop (Nr:1)  C--     end of thermodynamic k loop (Nr:1)
597          ENDDO          ENDDO
598    
599    C       All explicit advection/diffusion/sources should now be
600    C       done. The updated tracer field is in gPtr. Accumalate
601    C       explicit tendency and also reset gPtr to initial tracer
602    C       field for implicit matrix calculation
603    
604    #ifdef ALLOW_MATRIX
605            IF (useMATRIX)
606         &    CALL MATRIX_STORE_TENDENCY_EXP(bi,bj, myTime,myIter,myThid)
607    #endif
608    
609            iMin = 1
610            iMax = sNx
611            jMin = 1
612            jMax = sNy
613    
614  C--     Implicit vertical advection & diffusion  C--     Implicit vertical advection & diffusion
615    #ifndef ALLOW_OFFLINE
616            IF ( tempStepping .AND. implicitDiffusion ) THEN
617              CALL CALC_3D_DIFFUSIVITY(
618         I         bi,bj,iMin,iMax,jMin,jMax,
619         I         GAD_TEMPERATURE, useGMredi, useKPP,
620         O         kappaRk,
621         I         myThid)
622            ENDIF
623  #ifdef INCLUDE_IMPLVERTADV_CODE  #ifdef INCLUDE_IMPLVERTADV_CODE
624          IF ( tempImplVertAdv ) THEN          IF ( tempImplVertAdv ) THEN
625            CALL GAD_IMPLICIT_R(            CALL GAD_IMPLICIT_R(
626       I         tempImplVertAdv, tempAdvScheme, GAD_TEMPERATURE,       I         tempImplVertAdv, tempAdvScheme, GAD_TEMPERATURE,
627       I         kappaRT, wVel, theta,       I         kappaRk, wVel, theta,
628       U         gT,       U         gT,
629       I         bi, bj, myTime, myIter, myThid )       I         bi, bj, myTime, myIter, myThid )
630          ELSEIF ( tempStepping .AND. implicitDiffusion ) THEN          ELSEIF ( tempStepping .AND. implicitDiffusion ) THEN
# Line 789  C--     Implicit vertical advection & di Line 632  C--     Implicit vertical advection & di
632          IF     ( tempStepping .AND. implicitDiffusion ) THEN          IF     ( tempStepping .AND. implicitDiffusion ) THEN
633  #endif /* INCLUDE_IMPLVERTADV_CODE */  #endif /* INCLUDE_IMPLVERTADV_CODE */
634  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
635  CADJ STORE KappaRT(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte  CADJ STORE kappaRk(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
636  CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte  CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
637  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
638            CALL IMPLDIFF(            CALL IMPLDIFF(
639       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
640       I         deltaTtracer, KappaRT, recip_HFacC,       I         GAD_TEMPERATURE, kappaRk, recip_hFacC,
641       U         gT,       U         gT,
642       I         myThid )       I         myThid )
643          ENDIF          ENDIF
644    #endif /* ndef ALLOW_OFFLINE */
645    
646    #ifdef ALLOW_TIMEAVE
647            useVariableK = useKPP .OR. usePP81 .OR. useMY82 .OR. useGGL90
648         &       .OR. useGMredi .OR. ivdc_kappa.NE.0.
649            IF (taveFreq.GT.0. .AND. useVariableK ) THEN
650             IF (implicitDiffusion) THEN
651               CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, gT, kappaRk,
652         I                        Nr, 3, deltaTclock, bi, bj, myThid)
653    c        ELSE
654    c          CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, theta, kappaRT,
655    c    I                        Nr, 3, deltaTclock, bi, bj, myThid)
656             ENDIF
657            ENDIF
658    #endif /* ALLOW_TIMEAVE */
659    
660    #ifndef ALLOW_OFFLINE
661            IF ( saltStepping .AND. implicitDiffusion ) THEN
662              CALL CALC_3D_DIFFUSIVITY(
663         I         bi,bj,iMin,iMax,jMin,jMax,
664         I         GAD_SALINITY, useGMredi, useKPP,
665         O         kappaRk,
666         I         myThid)
667            ENDIF
668    
669  #ifdef INCLUDE_IMPLVERTADV_CODE  #ifdef INCLUDE_IMPLVERTADV_CODE
670          IF ( saltImplVertAdv ) THEN          IF ( saltImplVertAdv ) THEN
671            CALL GAD_IMPLICIT_R(            CALL GAD_IMPLICIT_R(
672       I         saltImplVertAdv, saltAdvScheme, GAD_SALINITY,       I         saltImplVertAdv, saltAdvScheme, GAD_SALINITY,
673       I         kappaRS, wVel, salt,       I         kappaRk, wVel, salt,
674       U         gS,       U         gS,
675       I         bi, bj, myTime, myIter, myThid )       I         bi, bj, myTime, myIter, myThid )
676          ELSEIF ( saltStepping .AND. implicitDiffusion ) THEN          ELSEIF ( saltStepping .AND. implicitDiffusion ) THEN
# Line 811  CADJ STORE gT(:,:,:,bi,bj) = comlev1_bib Line 678  CADJ STORE gT(:,:,:,bi,bj) = comlev1_bib
678          IF     ( saltStepping .AND. implicitDiffusion ) THEN          IF     ( saltStepping .AND. implicitDiffusion ) THEN
679  #endif /* INCLUDE_IMPLVERTADV_CODE */  #endif /* INCLUDE_IMPLVERTADV_CODE */
680  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
681  CADJ STORE KappaRS(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte  CADJ STORE kappaRk(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
682  CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte  CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
683  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
684            CALL IMPLDIFF(            CALL IMPLDIFF(
685       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
686       I         deltaTtracer, KappaRS, recip_HFacC,       I         GAD_SALINITY, kappaRk, recip_hFacC,
687       U         gS,       U         gS,
688       I         myThid )       I         myThid )
689          ENDIF          ENDIF
   
 #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  
690  #endif  #endif
691    
692  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
693  c #ifdef INCLUDE_IMPLVERTADV_CODE          IF     ( usePTRACERS ) THEN
694  c       IF ( usePTRACERS .AND. ptracerImplVertAdv ) THEN  C--     Vertical advection/diffusion (implicit) for passive tracers
695  c       ELSEIF ( usePTRACERS .AND. implicitDiffusion ) THEN             CALL PTRACERS_IMPLICIT(
696  c #else       U                             kappaRk,
697          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 )  
698          ENDIF          ENDIF
699  #endif /* ALLOW_PTRACERS */  #endif /* ALLOW_PTRACERS */
700    
# Line 857  C--      Apply open boundary conditions Line 710  C--      Apply open boundary conditions
710          ENDIF          ENDIF
711  #endif   /* ALLOW_OBCS */  #endif   /* ALLOW_OBCS */
712    
 #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 */  
   
713  #endif /* SINGLE_LAYER_MODE */  #endif /* SINGLE_LAYER_MODE */
714    
715  C--   end bi,bj loops.  C--   end bi,bj loops.
# Line 889  C--   end bi,bj loops. Line 723  C--   end bi,bj loops.
723         CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (THERMODYNAMICS)',myThid)
724         CALL DEBUG_STATS_RL(Nr,theta,'Theta (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,theta,'Theta (THERMODYNAMICS)',myThid)
725         CALL DEBUG_STATS_RL(Nr,salt,'Salt (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,salt,'Salt (THERMODYNAMICS)',myThid)
726         CALL DEBUG_STATS_RL(Nr,Gt,'Gt (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,gT,'Gt (THERMODYNAMICS)',myThid)
727         CALL DEBUG_STATS_RL(Nr,Gs,'Gs (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,gS,'Gs (THERMODYNAMICS)',myThid)
728         CALL DEBUG_STATS_RL(Nr,GtNm1,'GtNm1 (THERMODYNAMICS)',myThid)  #ifndef ALLOW_ADAMSBASHFORTH_3
729         CALL DEBUG_STATS_RL(Nr,GsNm1,'GsNm1 (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,gtNm1,'GtNm1 (THERMODYNAMICS)',myThid)
730           CALL DEBUG_STATS_RL(Nr,gsNm1,'GsNm1 (THERMODYNAMICS)',myThid)
731    #endif
732  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
733         IF ( usePTRACERS ) THEN         IF ( usePTRACERS ) THEN
734           CALL PTRACERS_DEBUG(myThid)           CALL PTRACERS_DEBUG(myThid)
735         ENDIF         ENDIF
736  #endif /* ALLOW_PTRACERS */  #endif /* ALLOW_PTRACERS */
737        ENDIF        ENDIF
738  #endif  #endif /* ALLOW_DEBUG */
739    
740  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
741           IF ( debugLevel .GE. debLevB )           IF ( debugLevel .GE. debLevB )
742       &    CALL DEBUG_LEAVE('THERMODYNAMICS',myThid)       &    CALL DEBUG_LEAVE('THERMODYNAMICS',myThid)
743  #endif  #endif
744    
745    #endif /* ALLOW_GENERIC_ADVDIFF */
746    
747        RETURN        RETURN
748        END        END
Legend:
Removed from v.1.65  
changed lines
  Added in v.1.99
  ViewVC Help
Powered by ViewVC 1.1.22