/[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.70 by jmc, Thu Jul 1 21:45:11 2004 UTC revision 1.119 by heimbach, Tue Oct 23 16:55:46 2007 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"
6  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_GENERIC_ADVDIFF
7  # include "PTRACERS_OPTIONS.h"  # include "GAD_OPTIONS.h"
8  #endif  #endif
9    
10  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 22  C     !INTERFACE: Line 22  C     !INTERFACE:
22        SUBROUTINE THERMODYNAMICS(myTime, myIter, myThid)        SUBROUTINE THERMODYNAMICS(myTime, myIter, myThid)
23  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
24  C     *==========================================================*  C     *==========================================================*
25  C     | SUBROUTINE THERMODYNAMICS                                  C     | SUBROUTINE THERMODYNAMICS
26  C     | o Controlling routine for the prognostic part of the        C     | o Controlling routine for the prognostic part of the
27  C     |   thermo-dynamics.                                          C     |   thermo-dynamics.
28  C     *===========================================================  C     *===========================================================
29  C     | The algorithm...  C     | The algorithm...
30  C     |  C     |
# Line 78  C     == Global variables === Line 78  C     == Global variables ===
78  #include "SIZE.h"  #include "SIZE.h"
79  #include "EEPARAMS.h"  #include "EEPARAMS.h"
80  #include "PARAMS.h"  #include "PARAMS.h"
81    #include "RESTART.h"
82  #include "DYNVARS.h"  #include "DYNVARS.h"
83  #include "GRID.h"  #include "GRID.h"
84  #include "GAD.h"  #ifdef ALLOW_GENERIC_ADVDIFF
85  #ifdef ALLOW_PASSIVE_TRACER  # include "GAD.h"
86  #include "TR1.h"  # include "GAD_SOM_VARS.h"
87  #endif  #endif
88  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
89  #include "PTRACERS.h"  # include "PTRACERS_SIZE.h"
90    # include "PTRACERS.h"
91  #endif  #endif
92  #ifdef ALLOW_TIMEAVE  #ifdef ALLOW_TIMEAVE
93  #include "TIMEAVE_STATV.h"  # include "TIMEAVE_STATV.h"
94  #endif  #endif
95    
96  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
97  # include "tamc.h"  # include "tamc.h"
98  # include "tamc_keys.h"  # include "tamc_keys.h"
99  # include "FFIELDS.h"  # include "FFIELDS.h"
100    # include "SURFACE.h"
101  # include "EOS.h"  # include "EOS.h"
102  # ifdef ALLOW_KPP  # ifdef ALLOW_KPP
103  #  include "KPP.h"  #  include "KPP.h"
# Line 116  C     myThid - Thread number for this in Line 119  C     myThid - Thread number for this in
119        INTEGER myIter        INTEGER myIter
120        INTEGER myThid        INTEGER myThid
121    
122    #ifdef ALLOW_GENERIC_ADVDIFF
123  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
124  C     == Local variables  C     == Local variables
125  C     xA, yA                 - Per block temporaries holding face areas  C     xA, yA                 - Per block temporaries holding face areas
126  C     uTrans, vTrans, rTrans - Per block temporaries holding flow  C     uFld, vFld, wFld       - Local copy of velocity field (3 components)
127  C                              transport  C     uTrans, vTrans, rTrans - Per block temporaries holding flow transport
128  C                              o uTrans: Zonal transport  C                              o uTrans: Zonal transport
129  C                              o vTrans: Meridional transport  C                              o vTrans: Meridional transport
130  C                              o rTrans: Vertical transport  C                              o rTrans: Vertical transport
# Line 130  C     fVer[STUV]               o fVer: V Line 134  C     fVer[STUV]               o fVer: V
134  C                                      is "pipelined" in the vertical  C                                      is "pipelined" in the vertical
135  C                                      so we need an fVer for each  C                                      so we need an fVer for each
136  C                                      variable.  C                                      variable.
137  C     rhoK, rhoKM1   - Density at current level, and level above  C     kappaRT,       - Total diffusion in vertical at level k, for T and S
138  C     KappaRT,       - Total diffusion in vertical for T and S.  C     kappaRS          (background + spatially varying, isopycnal term).
139  C     KappaRS          (background + spatially varying, isopycnal term).  C     kappaRTr       - Total diffusion in vertical at level k,
140    C                      for each passive Tracer
141    C     kappaRk        - Total diffusion in vertical, all levels, 1 tracer
142  C     useVariableK   = T when vertical diffusion is not constant  C     useVariableK   = T when vertical diffusion is not constant
143  C     iMin, iMax     - Ranges and sub-block indices on which calculations  C     iMin, iMax     - Ranges and sub-block indices on which calculations
144  C     jMin, jMax       are applied.  C     jMin, jMax       are applied.
# Line 142  C     kDown, km1       are switched with Line 148  C     kDown, km1       are switched with
148  C                      index into fVerTerm.  C                      index into fVerTerm.
149        _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
150        _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
151          _RL uFld    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
152          _RL vFld    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
153          _RL wFld    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
154        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
155        _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
156        _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 149  C                      index into fVerTe Line 158  C                      index into fVerTe
158        _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
159        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
160        _RL fVerS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
161  #ifdef ALLOW_PASSIVE_TRACER        _RL kappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
162        _RL fVerTr1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL kappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
 #endif  
163  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
164        _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)
165          _RL kappaRTr(1-Olx:sNx+Olx,1-Oly:sNy+Oly,PTRACERS_num)
166  #endif  #endif
167        _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 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)  
       _RL kp1Msk  
       LOGICAL useVariableK  
168        INTEGER iMin, iMax        INTEGER iMin, iMax
169        INTEGER jMin, jMax        INTEGER jMin, jMax
170        INTEGER bi, bj        INTEGER bi, bj
171        INTEGER i, j        INTEGER i, j
172        INTEGER k, km1, kup, kDown        INTEGER k, km1, kup, kDown
173    #ifdef ALLOW_ADAMSBASHFORTH_3
174          INTEGER iterNb, m1, m2
175    #endif
176    #ifdef ALLOW_TIMEAVE
177          LOGICAL useVariableK
178    #endif
179    #ifdef ALLOW_PTRACERS
180        INTEGER iTracer, ip        INTEGER iTracer, ip
181    #endif
182    
183  CEOP  CEOP
184    
185  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
186           IF ( debugLevel .GE. debLevB )           IF ( debugLevel .GE. debLevB )
187       &    CALL DEBUG_ENTER('THERMODYNAMICS',myThid)       &    CALL DEBUG_ENTER('THERMODYNAMICS',myThid)
188  #endif  #endif
189    
190  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
191  C--   dummy statement to end declaration part  C--   dummy statement to end declaration part
192        ikey = 1        ikey = 1
# Line 189  C--   HPF directive to help TAMC Line 198  C--   HPF directive to help TAMC
198  CHPF$ INDEPENDENT  CHPF$ INDEPENDENT
199  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
200    
201    C-- Compute correction at the surface for Lin Free Surf.
202    #ifdef ALLOW_AUTODIFF_TAMC
203          TsurfCor = 0. _d 0
204          SsurfCor = 0. _d 0
205    CADJ STORE theta,salt,wvel = comlev1, key = ikey_dynamics
206    #endif
207          IF (linFSConserveTr) THEN
208    #ifdef ALLOW_AUTODIFF_TAMC
209    CADJ STORE theta,salt,wvel = comlev1, key = ikey_dynamics
210    #endif
211           CALL CALC_WSURF_TR(theta,salt,wVel,
212         &                    myTime,myIter,myThid)
213          ENDIF
214    
215        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
216    
217  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
218  C--    HPF directive to help TAMC  C--    HPF directive to help TAMC
219  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS
220  CHPF$&                  ,utrans,vtrans,xA,yA  CHPF$&                  ,utrans,vtrans,xA,yA
221  CHPF$&                  ,KappaRT,KappaRS  CHPF$&                  ,kappaRT,kappaRS
222  CHPF$&                  )  CHPF$&                  )
223    # ifdef ALLOW_PTRACERS
224    CHPF$  INDEPENDENT, NEW (fVerP,kappaRTr)
225    # endif
226  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
227    
228         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
# Line 226  C     uninitialised but inert locations. Line 252  C     uninitialised but inert locations.
252            yA(i,j)        = 0. _d 0            yA(i,j)        = 0. _d 0
253            uTrans(i,j)    = 0. _d 0            uTrans(i,j)    = 0. _d 0
254            vTrans(i,j)    = 0. _d 0            vTrans(i,j)    = 0. _d 0
           rhok   (i,j)   = 0. _d 0  
           rhoKM1 (i,j)   = 0. _d 0  
255            rTrans (i,j)   = 0. _d 0            rTrans (i,j)   = 0. _d 0
256            rTransKp1(i,j) = 0. _d 0            rTransKp1(i,j) = 0. _d 0
257            fVerT  (i,j,1) = 0. _d 0            fVerT  (i,j,1) = 0. _d 0
258            fVerT  (i,j,2) = 0. _d 0            fVerT  (i,j,2) = 0. _d 0
259            fVerS  (i,j,1) = 0. _d 0            fVerS  (i,j,1) = 0. _d 0
260            fVerS  (i,j,2) = 0. _d 0            fVerS  (i,j,2) = 0. _d 0
261  #ifdef ALLOW_PASSIVE_TRACER            kappaRT(i,j)   = 0. _d 0
262            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  
263           ENDDO           ENDDO
264          ENDDO          ENDDO
265    
# Line 251  C     uninitialised but inert locations. Line 267  C     uninitialised but inert locations.
267           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
268            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
269  C This is currently also used by IVDC and Diagnostics  C This is currently also used by IVDC and Diagnostics
270             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  
            KappaRT(i,j,k)    = 0. _d 0  
            KappaRS(i,j,k)    = 0. _d 0  
271  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):
272             gT(i,j,k,bi,bj)   = 0. _d 0             gT(i,j,k,bi,bj)   = 0. _d 0
273             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 */  
274            ENDDO            ENDDO
275           ENDDO           ENDDO
276          ENDDO          ENDDO
277    
278          iMin = 1-OLx  #ifdef ALLOW_PTRACERS
279          iMax = sNx+OLx          IF ( usePTRACERS ) THEN
280          jMin = 1-OLy           DO ip=1,PTRACERS_num
281          jMax = sNy+OLy             DO j=1-OLy,sNy+OLy
282                DO i=1-OLx,sNx+OLx
283  #ifdef ALLOW_AUTODIFF_TAMC               fVerP  (i,j,1,ip) = 0. _d 0
284  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte               fVerP  (i,j,2,ip) = 0. _d 0
285  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte               kappaRTr(i,j,ip)  = 0. _d 0
286  CADJ STORE totphihyd              ENDDO
287  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte             ENDDO
288  #ifdef ALLOW_KPP           ENDDO
289  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  C-      set tracer tendency to zero:
290  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte           DO iTracer=1,PTRACERS_num
291  #endif            DO k=1,Nr
292  #endif /* ALLOW_AUTODIFF_TAMC */             DO j=1-OLy,sNy+OLy
293                DO i=1-OLx,sNx+OLx
294  #ifdef ALLOW_DEBUG               gPTr(i,j,k,bi,bj,itracer) = 0. _d 0
295          IF ( debugLevel .GE. debLevB )              ENDDO
296       &    CALL DEBUG_MSG('ENTERING UPWARD K LOOP',myThid)             ENDDO
297  #endif            ENDDO
298             ENDDO
299  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)  
300  #endif  #endif
             CALL CALC_IVDC(  
      I        bi, bj, iMin, iMax, jMin, jMax, k,  
      I        rhoKm1, rhoK,  
      I        myTime, myIter, myThid)  
           ENDIF  
301    
302  #endif /* SINGLE_LAYER_MODE */  #ifdef ALLOW_ADAMSBASHFORTH_3
303    C-      Apply AB on T,S :
304  C--     end of diagnostic k loop (Nr:1)          iterNb = myIter
305          ENDDO          IF (staggerTimeStep) iterNb = myIter - 1
306            m1 = 1 + MOD(iterNb+1,2)
307            m2 = 1 + MOD( iterNb ,2)
308    C       compute T^n+1/2 (stored in gtNm) extrapolating T forward in time
309            IF ( AdamsBashforth_T ) CALL ADAMS_BASHFORTH3(
310         I                                  bi, bj, 0,
311         U                                  theta, gtNm,
312         I                                  tempStartAB, iterNb, myThid )
313    C       compute S^n+1/2 (stored in gsNm) extrapolating S forward in time
314            IF ( AdamsBashforth_S ) CALL ADAMS_BASHFORTH3(
315         I                                  bi, bj, 0,
316         U                                  salt, gsNm,
317         I                                  saltStartAB, iterNb, myThid )
318    #endif /* ALLOW_ADAMSBASHFORTH_3 */
319    
320    c       iMin = 1-OLx
321    c       iMax = sNx+OLx
322    c       jMin = 1-OLy
323    c       jMax = sNy+OLy
324    
325  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
326  cph avoids recomputation of integrate_for_w  cph avoids recomputation of integrate_for_w
327  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
328  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
329    
 #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 */  
   
330  C--     Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h  C--     Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h
331  C--     MOST of THERMODYNAMICS will be disabled  C--     MOST of THERMODYNAMICS will be disabled
332  #ifndef SINGLE_LAYER_MODE  #ifndef SINGLE_LAYER_MODE
333    
 #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 */  
   
334  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
335  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
336  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
337  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
338  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
339  #ifdef ALLOW_PASSIVE_TRACER  # if ((defined ALLOW_DEPTH_CONTROL) || (defined NONLIN_FRSURF))
340  CADJ STORE tr1  (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE gtnm1(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
341  #endif  CADJ STORE gsnm1(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
342  #ifdef ALLOW_PTRACERS  # endif
 cph-- moved to forward_step to avoid key computation  
 cphCADJ STORE ptracer(:,:,:,bi,bj,itracer) = comlev1_bibj,  
 cphCADJ &                              key=itdkey, byte=isbyte  
 #endif  
343  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
344    
 #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 */  
   
345  #ifndef DISABLE_MULTIDIM_ADVECTION  #ifndef DISABLE_MULTIDIM_ADVECTION
346  C--     Some advection schemes are better calculated using a multi-dimensional  C--     Some advection schemes are better calculated using a multi-dimensional
347  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 567  C to be able to exclude this scheme to a Line 353  C to be able to exclude this scheme to a
353  C recomputation. It *is* differentiable, if you need it.  C recomputation. It *is* differentiable, if you need it.
354  C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to  C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to
355  C disable this section of code.  C disable this section of code.
356    #ifdef GAD_ALLOW_SOM_ADVECT
357            IF ( tempSOM_Advection ) THEN
358    #ifdef ALLOW_DEBUG
359              IF ( debugLevel .GE. debLevB )
360         &     CALL DEBUG_CALL('GAD_SOM_ADVECT',myThid)
361    #endif
362              CALL GAD_SOM_ADVECT(
363         I             tempImplVertAdv, tempAdvScheme, tempVertAdvScheme,
364         I             GAD_TEMPERATURE,
365         I             uVel, vVel, wVel, theta,
366         U             som_T,
367         O             gT,
368         I             bi,bj,myTime,myIter,myThid)
369            ELSEIF (tempMultiDimAdvec) THEN
370    #else /* GAD_ALLOW_SOM_ADVECT */
371          IF (tempMultiDimAdvec) THEN          IF (tempMultiDimAdvec) THEN
372    #endif /* GAD_ALLOW_SOM_ADVECT */
373  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
374            IF ( debugLevel .GE. debLevB )            IF ( debugLevel .GE. debLevB )
375       &     CALL DEBUG_CALL('GAD_ADVECTION',myThid)       &     CALL DEBUG_CALL('GAD_ADVECTION',myThid)
376  #endif  #endif
377            CALL GAD_ADVECTION(            CALL GAD_ADVECTION(
# Line 579  C disable this section of code. Line 381  C disable this section of code.
381       O             gT,       O             gT,
382       I             bi,bj,myTime,myIter,myThid)       I             bi,bj,myTime,myIter,myThid)
383          ENDIF          ENDIF
384    #ifdef GAD_ALLOW_SOM_ADVECT
385            IF ( saltSOM_Advection ) THEN
386    #ifdef ALLOW_DEBUG
387              IF ( debugLevel .GE. debLevB )
388         &     CALL DEBUG_CALL('GAD_SOM_ADVECT',myThid)
389    #endif
390              CALL GAD_SOM_ADVECT(
391         I             saltImplVertAdv, saltAdvScheme, saltVertAdvScheme,
392         I             GAD_SALINITY,
393         I             uVel, vVel, wVel, salt,
394         U             som_S,
395         O             gS,
396         I             bi,bj,myTime,myIter,myThid)
397            ELSEIF (saltMultiDimAdvec) THEN
398    #else /* GAD_ALLOW_SOM_ADVECT */
399          IF (saltMultiDimAdvec) THEN          IF (saltMultiDimAdvec) THEN
400    #endif /* GAD_ALLOW_SOM_ADVECT */
401  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
402            IF ( debugLevel .GE. debLevB )            IF ( debugLevel .GE. debLevB )
403       &     CALL DEBUG_CALL('GAD_ADVECTION',myThid)       &     CALL DEBUG_CALL('GAD_ADVECTION',myThid)
404  #endif  #endif
405            CALL GAD_ADVECTION(            CALL GAD_ADVECTION(
# Line 591  C disable this section of code. Line 409  C disable this section of code.
409       O             gS,       O             gS,
410       I             bi,bj,myTime,myIter,myThid)       I             bi,bj,myTime,myIter,myThid)
411          ENDIF          ENDIF
412    
413  C Since passive tracers are configurable separately from T,S we  C Since passive tracers are configurable separately from T,S we
414  C call the multi-dimensional method for PTRACERS regardless  C call the multi-dimensional method for PTRACERS regardless
415  C of whether multiDimAdvection is set or not.  C of whether multiDimAdvection is set or not.
416  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
417          IF ( usePTRACERS ) THEN          IF ( usePTRACERS ) THEN
418  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
419            IF ( debugLevel .GE. debLevB )            IF ( debugLevel .GE. debLevB )
420       &     CALL DEBUG_CALL('PTRACERS_ADVECTION',myThid)       &     CALL DEBUG_CALL('PTRACERS_ADVECTION',myThid)
421  #endif  #endif
422           CALL PTRACERS_ADVECTION( bi,bj,myIter,myTime,myThid )           CALL PTRACERS_ADVECTION( bi,bj,myIter,myTime,myThid )
# Line 606  C of whether multiDimAdvection is set or Line 425  C of whether multiDimAdvection is set or
425  #endif /* DISABLE_MULTIDIM_ADVECTION */  #endif /* DISABLE_MULTIDIM_ADVECTION */
426    
427  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
428          IF ( debugLevel .GE. debLevB )          IF ( debugLevel .GE. debLevB )
429       &    CALL DEBUG_MSG('ENTERING DOWNWARD K LOOP',myThid)       &    CALL DEBUG_MSG('ENTERING DOWNWARD K LOOP',myThid)
430  #endif  #endif
431    
# Line 615  C--     Start of thermodynamics loop Line 434  C--     Start of thermodynamics loop
434  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
435  C? Patrick Is this formula correct?  C? Patrick Is this formula correct?
436  cph Yes, but I rewrote it.  cph Yes, but I rewrote it.
437  cph Also, the KappaR? need the index and subscript k!  cph Also, the kappaR? need the index and subscript k!
438           kkey = (itdkey-1)*Nr + k           kkey = (itdkey-1)*Nr + k
439  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
440    
# Line 632  C--       kDown  Cycles through 2,1 to p Line 451  C--       kDown  Cycles through 2,1 to p
451            jMin = 1-OLy            jMin = 1-OLy
452            jMax = sNy+OLy            jMax = sNy+OLy
453    
454            kp1Msk=1.            IF (k.EQ.Nr) THEN
455            IF (k.EQ.Nr) kp1Msk=0.             DO j=1-Oly,sNy+Oly
456                DO i=1-Olx,sNx+Olx
457                 rTransKp1(i,j) = 0. _d 0
458                ENDDO
459               ENDDO
460              ELSE
461             DO j=1-Oly,sNy+Oly             DO j=1-Oly,sNy+Oly
462              DO i=1-Olx,sNx+Olx              DO i=1-Olx,sNx+Olx
463               rTransKp1(i,j) = kp1Msk*rTrans(i,j)               rTransKp1(i,j) = rTrans(i,j)
464              ENDDO              ENDDO
465             ENDDO             ENDDO
466              ENDIF
467  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
468  CADJ STORE rTransKp1(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE rTransKp1(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
469  #endif  #endif
470    
471  C--       Get temporary terms used by tendency routines  C--       Get temporary terms used by tendency routines :
472    C-        Calculate horizontal "volume transport" through tracer cell face
473    C         anelastic: uTrans,vTrans are scaled by rhoFacC (~ mass transport)
474            CALL CALC_COMMON_FACTORS (            CALL CALC_COMMON_FACTORS (
475       I         bi,bj,iMin,iMax,jMin,jMax,k,       I         uVel, vVel,
476       O         xA,yA,uTrans,vTrans,rTrans,maskUp,       O         uFld, vFld, uTrans, vTrans, xA, yA,
477       I         myThid)       I         k,bi,bj, myThid )
478    
479    C-        Calculate vertical "volume transport" through tracer cell face
480            IF (k.EQ.1) THEN            IF (k.EQ.1) THEN
481  C- Surface interface :  C-         Surface interface :
482             DO j=1-Oly,sNy+Oly             DO j=1-Oly,sNy+Oly
483              DO i=1-Olx,sNx+Olx              DO i=1-Olx,sNx+Olx
484               rTrans(i,j) = 0.               wFld(i,j)   = 0. _d 0
485                 maskUp(i,j) = 0. _d 0
486                 rTrans(i,j) = 0. _d 0
487              ENDDO              ENDDO
488             ENDDO             ENDDO
489            ELSE            ELSE
490  C- Interior interface :  C-         Interior interface :
491    C          anelastic: rTrans is scaled by rhoFacF (~ mass transport)
492             DO j=1-Oly,sNy+Oly             DO j=1-Oly,sNy+Oly
493              DO i=1-Olx,sNx+Olx              DO i=1-Olx,sNx+Olx
494               rTrans(i,j) = rTrans(i,j)*maskC(i,j,k-1,bi,bj)               wFld(i,j)   = wVel(i,j,k,bi,bj)
495                 maskUp(i,j) = maskC(i,j,k-1,bi,bj)*maskC(i,j,k,bi,bj)
496                 rTrans(i,j) = wFld(i,j)*rA(i,j,bi,bj)*maskUp(i,j)
497         &                              *deepFac2F(k)*rhoFacF(k)
498              ENDDO              ENDDO
499             ENDDO             ENDDO
500            ENDIF            ENDIF
501    
502  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
   
503  C--   Residual transp = Bolus transp + Eulerian transp  C--   Residual transp = Bolus transp + Eulerian transp
504            IF (useGMRedi) THEN            IF (useGMRedi) THEN
505              CALL GMREDI_CALC_UVFLOW(              CALL GMREDI_CALC_UVFLOW(
506       &            uTrans, vTrans, bi, bj, k, myThid)       U                  uFld, vFld, uTrans, vTrans,
507              IF (K.GE.2) CALL GMREDI_CALC_WFLOW(       I                  k, bi, bj, myThid )
508       &                    rTrans, bi, bj, k, myThid)              IF (K.GE.2) THEN
509                  CALL GMREDI_CALC_WFLOW(
510         U                  wFld, rTrans,
511         I                  k, bi, bj, myThid )
512                ENDIF
513            ENDIF            ENDIF
514    # ifdef ALLOW_AUTODIFF_TAMC
 #ifdef ALLOW_AUTODIFF_TAMC  
515  CADJ STORE rTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE rTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
516  #ifdef GM_BOLUS_ADVEC  CADJ STORE wfld(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
517    # ifdef GM_BOLUS_ADVEC
518    CADJ STORE ufld(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
519    CADJ STORE vfld(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
520  CADJ STORE uTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE uTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
521  CADJ STORE vTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE vTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
522  #endif  # endif
523  #endif /* ALLOW_AUTODIFF_TAMC */  # endif /* ALLOW_AUTODIFF_TAMC */
   
524  #endif /* ALLOW_GMREDI */  #endif /* ALLOW_GMREDI */
525    
526  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL
527  C--      Calculate the total vertical diffusivity  C--      Calculate the total vertical diffusivity
528           CALL CALC_DIFFUSIVITY(            IF ( .NOT.implicitDiffusion ) THEN
529       I        bi,bj,iMin,iMax,jMin,jMax,k,              CALL CALC_DIFFUSIVITY(
530       I        maskUp,       I          bi,bj,iMin,iMax,jMin,jMax,k,
531       O        KappaRT,KappaRS,       I          maskUp,
532       I        myThid)       O          kappaRT,kappaRS,
533         I          myThid)
534              ENDIF
535  # ifdef ALLOW_AUTODIFF_TAMC  # ifdef ALLOW_AUTODIFF_TAMC
536  CADJ STORE KappaRT(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE kappaRT(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
537  CADJ STORE KappaRS(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE kappaRS(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
538  # endif /* ALLOW_AUTODIFF_TAMC */  # endif /* ALLOW_AUTODIFF_TAMC */
539  #endif  #endif
540    
# Line 704  CADJ STORE KappaRS(:,:,k)    = comlev1_b Line 544  CADJ STORE KappaRS(:,:,k)    = comlev1_b
544            jMax = sNy+OLy-1            jMax = sNy+OLy-1
545    
546  C--      Calculate active tracer tendencies (gT,gS,...)  C--      Calculate active tracer tendencies (gT,gS,...)
547  C        and step forward storing result in gTnm1, gSnm1, etc.  C        and step forward storing result in gT, gS, etc.
548    C--
549    # ifdef ALLOW_AUTODIFF_TAMC
550    #  if ((defined NONLIN_FRSURF) || (defined ALLOW_DEPTH_CONTROL)) && (defined ALLOW_GMREDI)
551    #   ifdef GM_NON_UNITY_DIAGONAL
552    CADJ STORE kux(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
553    CADJ STORE kvy(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
554    #   endif
555    #   ifdef GM_EXTRA_DIAGONAL
556    CADJ STORE kuz(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
557    CADJ STORE kvz(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
558    #   endif
559    #  endif
560    # endif /* ALLOW_AUTODIFF_TAMC */
561    C
562    #ifdef ALLOW_AUTODIFF_TAMC
563    # if (defined NONLIN_FRSURF) || (defined ALLOW_DEPTH_CONTROL)
564    cph-test
565    CADJ STORE uFld(:,:), vFld(:,:), wFld(:,:)
566    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
567    CADJ STORE uTrans(:,:), vTrans(:,:)
568    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
569    CADJ STORE xA(:,:), yA(:,:)
570    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
571    # endif
572    #endif /* ALLOW_AUTODIFF_TAMC */
573    C
574           IF ( tempStepping ) THEN           IF ( tempStepping ) THEN
575             CALL CALC_GT(  #ifdef ALLOW_AUTODIFF_TAMC
576    CADJ STORE gTnm1(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
577    # if (defined NONLIN_FRSURF) || (defined ALLOW_DEPTH_CONTROL)
578    CADJ STORE gt(:,:,:,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
579    CADJ STORE fvert(:,:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
580    # endif
581    #endif /* ALLOW_AUTODIFF_TAMC */
582              CALL CALC_GT(
583       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
584       I         xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,       I         xA, yA, maskUp, uFld, vFld, wFld,
585       I         KappaRT,       I         uTrans, vTrans, rTrans, rTransKp1,
586         I         kappaRT,
587       U         fVerT,       U         fVerT,
588       I         myTime,myIter,myThid)       I         myTime,myIter,myThid)
589    #ifdef ALLOW_ADAMSBASHFORTH_3
590              IF ( AdamsBashforth_T ) THEN
591             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
592       I         bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,       I         bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,
593       I         theta, gT,       I         gtNm(1-Olx,1-Oly,1,1,1,m2),
594         U         gT,
595       I         myIter, myThid)       I         myIter, myThid)
596              ELSE
597    #endif
598               CALL TIMESTEP_TRACER(
599         I         bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,
600         I         theta,
601         U         gT,
602         I         myIter, myThid)
603    #ifdef ALLOW_ADAMSBASHFORTH_3
604              ENDIF
605    #endif
606           ENDIF           ENDIF
607    
608           IF ( saltStepping ) THEN           IF ( saltStepping ) THEN
609             CALL CALC_GS(  #ifdef ALLOW_AUTODIFF_TAMC
610    CADJ STORE gSnm1(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
611    # if (defined NONLIN_FRSURF) || (defined ALLOW_DEPTH_CONTROL)
612    CADJ STORE gs(:,:,:,bi,bj)    = comlev1_bibj_k, key=kkey, byte=isbyte
613    CADJ STORE fvers(:,:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
614    # endif
615    #endif /* ALLOW_AUTODIFF_TAMC */
616              CALL CALC_GS(
617       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
618       I         xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,       I         xA, yA, maskUp, uFld, vFld, wFld,
619       I         KappaRS,       I         uTrans, vTrans, rTrans, rTransKp1,
620         I         kappaRS,
621       U         fVerS,       U         fVerS,
622       I         myTime,myIter,myThid)       I         myTime,myIter,myThid)
623    #ifdef ALLOW_ADAMSBASHFORTH_3
624              IF ( AdamsBashforth_S ) THEN
625             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
626       I         bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,       I         bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,
627       I         salt, gS,       I         gsNm(1-Olx,1-Oly,1,1,1,m2),
628         U         gS,
629       I         myIter, myThid)       I         myIter, myThid)
630           ENDIF            ELSE
631  #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)  
632             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
633       I         bi,bj,iMin,iMax,jMin,jMax,k,tracerAdvScheme,       I         bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,
634       I         Tr1, gTr1,       I         salt,
635       I         myIter,myThid)       U         gS,
636           ENDIF       I         myIter, myThid)
637    #ifdef ALLOW_ADAMSBASHFORTH_3
638              ENDIF
639  #endif  #endif
640             ENDIF
641    
642  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
643           IF ( usePTRACERS ) THEN           IF ( usePTRACERS ) THEN
644               IF ( .NOT.implicitDiffusion ) THEN
645                 CALL PTRACERS_CALC_DIFF(
646         I            bi,bj,iMin,iMax,jMin,jMax,k,
647         I            maskUp,
648         O            kappaRTr,
649         I            myThid)
650               ENDIF
651    # ifdef ALLOW_AUTODIFF_TAMC
652    CADJ STORE kappaRTr(:,:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
653    # endif /* ALLOW_AUTODIFF_TAMC */
654             CALL PTRACERS_INTEGRATE(             CALL PTRACERS_INTEGRATE(
655       I         bi,bj,k,       I          bi,bj,k,
656       I         xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,       I          xA, yA, maskUp, uFld, vFld, wFld,
657       X         fVerP, KappaRS,       I          uTrans, vTrans, rTrans, rTransKp1,
658       I         myIter,myTime,myThid)       I          kappaRTr,
659         U          fVerP,
660         I          myTime,myIter,myThid)
661           ENDIF           ENDIF
662  #endif /* ALLOW_PTRACERS */  #endif /* ALLOW_PTRACERS */
663    
# Line 777  CADJ &   , key = kkey, byte = isbyte Line 683  CADJ &   , key = kkey, byte = isbyte
683  C--     end of thermodynamic k loop (Nr:1)  C--     end of thermodynamic k loop (Nr:1)
684          ENDDO          ENDDO
685    
686    C       All explicit advection/diffusion/sources should now be
687    C       done. The updated tracer field is in gPtr. Accumalate
688    C       explicit tendency and also reset gPtr to initial tracer
689    C       field for implicit matrix calculation
690    
691    #ifdef ALLOW_MATRIX
692            IF (useMATRIX)
693         &    CALL MATRIX_STORE_TENDENCY_EXP(bi,bj, myTime,myIter,myThid)
694    #endif
695    
696            iMin = 1
697            iMax = sNx
698            jMin = 1
699            jMax = sNy
700    
701  C--     Implicit vertical advection & diffusion  C--     Implicit vertical advection & diffusion
702            IF ( tempStepping .AND. implicitDiffusion ) THEN
703              CALL CALC_3D_DIFFUSIVITY(
704         I         bi,bj,iMin,iMax,jMin,jMax,
705         I         GAD_TEMPERATURE, useGMredi, useKPP,
706         O         kappaRk,
707         I         myThid)
708            ENDIF
709  #ifdef INCLUDE_IMPLVERTADV_CODE  #ifdef INCLUDE_IMPLVERTADV_CODE
710          IF ( tempImplVertAdv ) THEN          IF ( tempImplVertAdv ) THEN
711            CALL GAD_IMPLICIT_R(            CALL GAD_IMPLICIT_R(
712       I         tempImplVertAdv, tempAdvScheme, GAD_TEMPERATURE,       I         tempImplVertAdv, tempAdvScheme, GAD_TEMPERATURE,
713       I         kappaRT, wVel, theta,       I         kappaRk, wVel, theta,
714       U         gT,       U         gT,
715       I         bi, bj, myTime, myIter, myThid )       I         bi, bj, myTime, myIter, myThid )
716          ELSEIF ( tempStepping .AND. implicitDiffusion ) THEN          ELSEIF ( tempStepping .AND. implicitDiffusion ) THEN
# Line 791  C--     Implicit vertical advection & di Line 718  C--     Implicit vertical advection & di
718          IF     ( tempStepping .AND. implicitDiffusion ) THEN          IF     ( tempStepping .AND. implicitDiffusion ) THEN
719  #endif /* INCLUDE_IMPLVERTADV_CODE */  #endif /* INCLUDE_IMPLVERTADV_CODE */
720  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
721  CADJ STORE KappaRT(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte  CADJ STORE kappaRk(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
722  CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte  CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
723  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
724            CALL IMPLDIFF(            CALL IMPLDIFF(
725       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
726       I         deltaTtracer, KappaRT, recip_HFacC,       I         GAD_TEMPERATURE, kappaRk, recip_hFacC,
727       U         gT,       U         gT,
728       I         myThid )       I         myThid )
729          ENDIF          ENDIF
730    
731    #ifdef ALLOW_TIMEAVE
732            useVariableK = useKPP .OR. usePP81 .OR. useMY82 .OR. useGGL90
733         &       .OR. useGMredi .OR. ivdc_kappa.NE.0.
734            IF (taveFreq.GT.0. .AND. useVariableK ) THEN
735             IF (implicitDiffusion) THEN
736               CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, gT, kappaRk,
737         I                        Nr, 3, deltaTclock, bi, bj, myThid)
738    c        ELSE
739    c          CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, theta, kappaRT,
740    c    I                        Nr, 3, deltaTclock, bi, bj, myThid)
741             ENDIF
742            ENDIF
743    #endif /* ALLOW_TIMEAVE */
744    
745            IF ( saltStepping .AND. implicitDiffusion ) THEN
746              CALL CALC_3D_DIFFUSIVITY(
747         I         bi,bj,iMin,iMax,jMin,jMax,
748         I         GAD_SALINITY, useGMredi, useKPP,
749         O         kappaRk,
750         I         myThid)
751            ENDIF
752    
753  #ifdef INCLUDE_IMPLVERTADV_CODE  #ifdef INCLUDE_IMPLVERTADV_CODE
754          IF ( saltImplVertAdv ) THEN          IF ( saltImplVertAdv ) THEN
755            CALL GAD_IMPLICIT_R(            CALL GAD_IMPLICIT_R(
756       I         saltImplVertAdv, saltAdvScheme, GAD_SALINITY,       I         saltImplVertAdv, saltAdvScheme, GAD_SALINITY,
757       I         kappaRS, wVel, salt,       I         kappaRk, wVel, salt,
758       U         gS,       U         gS,
759       I         bi, bj, myTime, myIter, myThid )       I         bi, bj, myTime, myIter, myThid )
760          ELSEIF ( saltStepping .AND. implicitDiffusion ) THEN          ELSEIF ( saltStepping .AND. implicitDiffusion ) THEN
# Line 813  CADJ STORE gT(:,:,:,bi,bj) = comlev1_bib Line 762  CADJ STORE gT(:,:,:,bi,bj) = comlev1_bib
762          IF     ( saltStepping .AND. implicitDiffusion ) THEN          IF     ( saltStepping .AND. implicitDiffusion ) THEN
763  #endif /* INCLUDE_IMPLVERTADV_CODE */  #endif /* INCLUDE_IMPLVERTADV_CODE */
764  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
765  CADJ STORE KappaRS(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte  CADJ STORE kappaRk(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
766  CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte  CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
767  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
768            CALL IMPLDIFF(            CALL IMPLDIFF(
769       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
770       I         deltaTtracer, KappaRS, recip_HFacC,       I         GAD_SALINITY, kappaRk, recip_hFacC,
771       U         gS,       U         gS,
772       I         myThid )       I         myThid )
773          ENDIF          ENDIF
774    
 #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  
   
775  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
776  c #ifdef INCLUDE_IMPLVERTADV_CODE          IF     ( usePTRACERS ) THEN
777  c       IF ( usePTRACERS .AND. ptracerImplVertAdv ) THEN  C--     Vertical advection/diffusion (implicit) for passive tracers
778  c       ELSEIF ( usePTRACERS .AND. implicitDiffusion ) THEN             CALL PTRACERS_IMPLICIT(
779  c #else       U                             kappaRk,
780          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 )  
781          ENDIF          ENDIF
782  #endif /* ALLOW_PTRACERS */  #endif /* ALLOW_PTRACERS */
783    
# Line 859  C--      Apply open boundary conditions Line 793  C--      Apply open boundary conditions
793          ENDIF          ENDIF
794  #endif   /* ALLOW_OBCS */  #endif   /* ALLOW_OBCS */
795    
 #ifdef ALLOW_TIMEAVE  
 #ifndef HRCUBE  
         IF (taveFreq.GT.0. .AND. ivdc_kappa.NE.0.) THEN  
           CALL TIMEAVE_CUMULATE(ConvectCountTave, IVDConvCount,  
      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 */  
   
796  #endif /* SINGLE_LAYER_MODE */  #endif /* SINGLE_LAYER_MODE */
797    
798  C--   end bi,bj loops.  C--   end bi,bj loops.
# Line 891  C--   end bi,bj loops. Line 806  C--   end bi,bj loops.
806         CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (THERMODYNAMICS)',myThid)
807         CALL DEBUG_STATS_RL(Nr,theta,'Theta (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,theta,'Theta (THERMODYNAMICS)',myThid)
808         CALL DEBUG_STATS_RL(Nr,salt,'Salt (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,salt,'Salt (THERMODYNAMICS)',myThid)
809         CALL DEBUG_STATS_RL(Nr,Gt,'Gt (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,gT,'Gt (THERMODYNAMICS)',myThid)
810         CALL DEBUG_STATS_RL(Nr,Gs,'Gs (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,gS,'Gs (THERMODYNAMICS)',myThid)
811         CALL DEBUG_STATS_RL(Nr,GtNm1,'GtNm1 (THERMODYNAMICS)',myThid)  #ifndef ALLOW_ADAMSBASHFORTH_3
812         CALL DEBUG_STATS_RL(Nr,GsNm1,'GsNm1 (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,gtNm1,'GtNm1 (THERMODYNAMICS)',myThid)
813           CALL DEBUG_STATS_RL(Nr,gsNm1,'GsNm1 (THERMODYNAMICS)',myThid)
814    #endif
815  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
816         IF ( usePTRACERS ) THEN         IF ( usePTRACERS ) THEN
817           CALL PTRACERS_DEBUG(myThid)           CALL PTRACERS_DEBUG(myThid)
818         ENDIF         ENDIF
819  #endif /* ALLOW_PTRACERS */  #endif /* ALLOW_PTRACERS */
820        ENDIF        ENDIF
821  #endif  #endif /* ALLOW_DEBUG */
822    
823  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
824           IF ( debugLevel .GE. debLevB )           IF ( debugLevel .GE. debLevB )
825       &    CALL DEBUG_LEAVE('THERMODYNAMICS',myThid)       &    CALL DEBUG_LEAVE('THERMODYNAMICS',myThid)
826  #endif  #endif
827    
828    #endif /* ALLOW_GENERIC_ADVDIFF */
829    
830        RETURN        RETURN
831        END        END
Legend:
Removed from v.1.70  
changed lines
  Added in v.1.119
  ViewVC Help
Powered by ViewVC 1.1.22