/[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.50 by jmc, Tue Oct 7 04:31:30 2003 UTC revision 1.109 by jmc, Tue Dec 5 05:25:08 2006 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    
4    #include "PACKAGES_CONFIG.h"
5  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
6    
7  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
8  # ifdef ALLOW_GMREDI  # ifdef ALLOW_GMREDI
9  #  include "GMREDI_OPTIONS.h"  #  include "GMREDI_OPTIONS.h"
# Line 9  C $Name$ Line 11  C $Name$
11  # ifdef ALLOW_KPP  # ifdef ALLOW_KPP
12  #  include "KPP_OPTIONS.h"  #  include "KPP_OPTIONS.h"
13  # endif  # endif
 #ifdef ALLOW_PTRACERS  
 # include "PTRACERS_OPTIONS.h"  
 #endif  
14  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
15    
16  CBOP  CBOP
# Line 20  C     !INTERFACE: Line 19  C     !INTERFACE:
19        SUBROUTINE THERMODYNAMICS(myTime, myIter, myThid)        SUBROUTINE THERMODYNAMICS(myTime, myIter, myThid)
20  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
21  C     *==========================================================*  C     *==========================================================*
22  C     | SUBROUTINE THERMODYNAMICS                                  C     | SUBROUTINE THERMODYNAMICS
23  C     | o Controlling routine for the prognostic part of the        C     | o Controlling routine for the prognostic part of the
24  C     |   thermo-dynamics.                                          C     |   thermo-dynamics.
25  C     *===========================================================  C     *===========================================================
26  C     | The algorithm...  C     | The algorithm...
27  C     |  C     |
# Line 78  C     == Global variables === Line 77  C     == Global variables ===
77  #include "PARAMS.h"  #include "PARAMS.h"
78  #include "DYNVARS.h"  #include "DYNVARS.h"
79  #include "GRID.h"  #include "GRID.h"
80    #ifdef ALLOW_GENERIC_ADVDIFF
81  #include "GAD.h"  #include "GAD.h"
 #ifdef ALLOW_PASSIVE_TRACER  
 #include "TR1.h"  
82  #endif  #endif
83  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
84    #include "PTRACERS_SIZE.h"
85  #include "PTRACERS.h"  #include "PTRACERS.h"
86  #endif  #endif
87  #ifdef ALLOW_TIMEAVE  #ifdef ALLOW_TIMEAVE
# Line 100  C     == Global variables === Line 99  C     == Global variables ===
99  # ifdef ALLOW_GMREDI  # ifdef ALLOW_GMREDI
100  #  include "GMREDI.h"  #  include "GMREDI.h"
101  # endif  # endif
102    # ifdef ALLOW_EBM
103    #  include "EBM.h"
104    # endif
105  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
106    
107  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
# Line 111  C     myThid - Thread number for this in Line 113  C     myThid - Thread number for this in
113        INTEGER myIter        INTEGER myIter
114        INTEGER myThid        INTEGER myThid
115    
116    #ifdef ALLOW_GENERIC_ADVDIFF
117  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
118  C     == Local variables  C     == Local variables
119  C     xA, yA                 - Per block temporaries holding face areas  C     xA, yA                 - Per block temporaries holding face areas
120  C     uTrans, vTrans, rTrans - Per block temporaries holding flow  C     uFld, vFld, wFld       - Local copy of velocity field (3 components)
121  C                              transport  C     uTrans, vTrans, rTrans - Per block temporaries holding flow transport
122  C                              o uTrans: Zonal transport  C                              o uTrans: Zonal transport
123  C                              o vTrans: Meridional transport  C                              o vTrans: Meridional transport
124  C                              o rTrans: Vertical transport  C                              o rTrans: Vertical transport
125    C     rTransKp1                o vertical volume transp. at interface k+1
126  C     maskUp                   o maskUp: land/water mask for W points  C     maskUp                   o maskUp: land/water mask for W points
127  C     fVer[STUV]               o fVer: Vertical flux term - note fVer  C     fVer[STUV]               o fVer: Vertical flux term - note fVer
128  C                                      is "pipelined" in the vertical  C                                      is "pipelined" in the vertical
129  C                                      so we need an fVer for each  C                                      so we need an fVer for each
130  C                                      variable.  C                                      variable.
131  C     rhoK, rhoKM1   - Density at current level, and level above  C     kappaRT,       - Total diffusion in vertical at level k, for T and S
132  C     phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean)  C     kappaRS          (background + spatially varying, isopycnal term).
133  C     phiSurfY             or geopotentiel (atmos) in X and Y direction  C     kappaRTr       - Total diffusion in vertical at level k,
134  C     KappaRT,       - Total diffusion in vertical for T and S.  C                      for each passive Tracer
135  C     KappaRS          (background + spatially varying, isopycnal term).  C     kappaRk        - Total diffusion in vertical, all levels, 1 tracer
136  C     useVariableK   = T when vertical diffusion is not constant  C     useVariableK   = T when vertical diffusion is not constant
137  C     iMin, iMax     - Ranges and sub-block indices on which calculations  C     iMin, iMax     - Ranges and sub-block indices on which calculations
138  C     jMin, jMax       are applied.  C     jMin, jMax       are applied.
# Line 138  C     kDown, km1       are switched with Line 142  C     kDown, km1       are switched with
142  C                      index into fVerTerm.  C                      index into fVerTerm.
143        _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
144        _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
145          _RL uFld    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
146          _RL vFld    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
147          _RL wFld    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
148        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
149        _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
150        _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
151          _RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
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        _RL fVerTr1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL kappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
156        _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL kappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
157        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  #ifdef ALLOW_PTRACERS
158        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fVerP   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2,PTRACERS_num)
159        _RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL kappaRTr(1-Olx:sNx+Olx,1-Oly:sNy+Oly,PTRACERS_num)
160        _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)  #endif
161        _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL kappaRk (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)  
       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        INTEGER iTracer  #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
175    #endif
176    
177  CEOP  CEOP
178    
179  #ifndef DISABLE_DEBUGMODE  #ifdef ALLOW_DEBUG
180           IF ( debugLevel .GE. debLevB )           IF ( debugLevel .GE. debLevB )
181       &    CALL DEBUG_ENTER('FORWARD_STEP',myThid)       &    CALL DEBUG_ENTER('THERMODYNAMICS',myThid)
182  #endif  #endif
183    
184  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
185  C--   dummy statement to end declaration part  C--   dummy statement to end declaration part
186        ikey = 1        ikey = 1
# Line 188  CHPF$ INDEPENDENT Line 198  CHPF$ INDEPENDENT
198  C--    HPF directive to help TAMC  C--    HPF directive to help TAMC
199  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS
200  CHPF$&                  ,utrans,vtrans,xA,yA  CHPF$&                  ,utrans,vtrans,xA,yA
201  CHPF$&                  ,KappaRT,KappaRS  CHPF$&                  ,kappaRT,kappaRS
202  CHPF$&                  )  CHPF$&                  )
203    # ifdef ALLOW_PTRACERS
204    CHPF$  INDEPENDENT, NEW (fVerP,kappaRTr)
205    # endif
206  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
207    
208         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
# Line 219  C     uninitialised but inert locations. Line 232  C     uninitialised but inert locations.
232            yA(i,j)        = 0. _d 0            yA(i,j)        = 0. _d 0
233            uTrans(i,j)    = 0. _d 0            uTrans(i,j)    = 0. _d 0
234            vTrans(i,j)    = 0. _d 0            vTrans(i,j)    = 0. _d 0
           rhok   (i,j)   = 0. _d 0  
           rhoKM1 (i,j)   = 0. _d 0  
           phiSurfX(i,j)  = 0. _d 0  
           phiSurfY(i,j)  = 0. _d 0  
235            rTrans (i,j)   = 0. _d 0            rTrans (i,j)   = 0. _d 0
236              rTransKp1(i,j) = 0. _d 0
237            fVerT  (i,j,1) = 0. _d 0            fVerT  (i,j,1) = 0. _d 0
238            fVerT  (i,j,2) = 0. _d 0            fVerT  (i,j,2) = 0. _d 0
239            fVerS  (i,j,1) = 0. _d 0            fVerS  (i,j,1) = 0. _d 0
240            fVerS  (i,j,2) = 0. _d 0            fVerS  (i,j,2) = 0. _d 0
241            fVerTr1(i,j,1) = 0. _d 0            kappaRT(i,j)   = 0. _d 0
242            fVerTr1(i,j,2) = 0. _d 0            kappaRS(i,j)   = 0. _d 0
243           ENDDO           ENDDO
244          ENDDO          ENDDO
245    
# Line 237  C     uninitialised but inert locations. Line 247  C     uninitialised but inert locations.
247           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
248            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
249  C This is currently also used by IVDC and Diagnostics  C This is currently also used by IVDC and Diagnostics
250             sigmaX(i,j,k) = 0. _d 0             kappaRk(i,j,k)    = 0. _d 0
            sigmaY(i,j,k) = 0. _d 0  
            sigmaR(i,j,k) = 0. _d 0  
            ConvectCount(i,j,k) = 0.  
            KappaRT(i,j,k)    = 0. _d 0  
            KappaRS(i,j,k)    = 0. _d 0  
251  C- tracer tendency needs to be set to zero (moved here from gad_calc_rhs):  C- tracer tendency needs to be set to zero (moved here from gad_calc_rhs):
252             gT(i,j,k,bi,bj)   = 0. _d 0             gT(i,j,k,bi,bj)   = 0. _d 0
253             gS(i,j,k,bi,bj)   = 0. _d 0             gS(i,j,k,bi,bj)   = 0. _d 0
 # ifdef ALLOW_PASSIVE_TRACER  
            gTr1(i,j,k,bi,bj) = 0. _d 0  
 # endif  
 # ifdef ALLOW_PTRACERS  
            DO iTracer=1,PTRACERS_numInUse  
             gPTr(i,j,k,bi,bj,itracer) = 0. _d 0  
            ENDDO  
 # endif  
 #ifdef ALLOW_AUTODIFF_TAMC  
 cph all the following init. are necessary for TAF  
 cph although some of these are re-initialised later.  
 # ifdef ALLOW_GMREDI  
            Kwx(i,j,k,bi,bj)  = 0. _d 0  
            Kwy(i,j,k,bi,bj)  = 0. _d 0  
            Kwz(i,j,k,bi,bj)  = 0. _d 0  
 #  ifdef GM_NON_UNITY_DIAGONAL  
            Kux(i,j,k,bi,bj)  = 0. _d 0  
            Kvy(i,j,k,bi,bj)  = 0. _d 0  
 #  endif  
 #  ifdef GM_EXTRA_DIAGONAL  
            Kuz(i,j,k,bi,bj)  = 0. _d 0  
            Kvz(i,j,k,bi,bj)  = 0. _d 0  
 #  endif  
 #  ifdef GM_BOLUS_ADVEC  
            GM_PsiX(i,j,k,bi,bj)  = 0. _d 0  
            GM_PsiY(i,j,k,bi,bj)  = 0. _d 0  
 #  endif  
 #  ifdef GM_VISBECK_VARIABLE_K  
            VisbeckK(i,j,bi,bj)   = 0. _d 0  
 #  endif  
 # endif /* ALLOW_GMREDI */  
 #endif /* ALLOW_AUTODIFF_TAMC */  
254            ENDDO            ENDDO
255           ENDDO           ENDDO
256          ENDDO          ENDDO
257    
258          iMin = 1-OLx  #ifdef ALLOW_PTRACERS
259          iMax = sNx+OLx          IF ( usePTRACERS ) THEN
260          jMin = 1-OLy           DO ip=1,PTRACERS_num
261          jMax = sNy+OLy             DO j=1-OLy,sNy+OLy
262                DO i=1-OLx,sNx+OLx
263  #ifdef ALLOW_AUTODIFF_TAMC               fVerP  (i,j,1,ip) = 0. _d 0
264  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte               fVerP  (i,j,2,ip) = 0. _d 0
265  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte               kappaRTr(i,j,ip)  = 0. _d 0
266  CADJ STORE totphihyd              ENDDO
267  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte             ENDDO
268  #ifdef ALLOW_KPP           ENDDO
269  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  C-      set tracer tendency to zero:
270  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte           DO iTracer=1,PTRACERS_num
271  #endif            DO k=1,Nr
272  #endif /* ALLOW_AUTODIFF_TAMC */             DO j=1-OLy,sNy+OLy
273                DO i=1-OLx,sNx+OLx
274  #ifndef DISABLE_DEBUGMODE               gPTr(i,j,k,bi,bj,itracer) = 0. _d 0
275          IF ( debugLevel .GE. debLevB )              ENDDO
276       &    CALL DEBUG_MSG('ENTERING UPWARD K LOOP',myThid)             ENDDO
277              ENDDO
278             ENDDO
279            ENDIF
280  #endif  #endif
281    
282  C--     Start of diagnostic loop  #ifdef ALLOW_ADAMSBASHFORTH_3
283          DO k=Nr,1,-1  C-      Apply AB on T,S :
284            iterNb = myIter
285  #ifdef ALLOW_AUTODIFF_TAMC          IF (staggerTimeStep) iterNb = myIter - 1
286  C? Patrick, is this formula correct now that we change the loop range?          m1 = 1 + MOD(iterNb+1,2)
287  C? Do we still need this?          m2 = 1 + MOD( iterNb ,2)
288  cph kkey formula corrected.  C       compute T^n+1/2 (stored in gtNm) extrapolating T forward in time
289  cph Needed for rhok, rhokm1, in the case useGMREDI.          IF ( AdamsBashforth_T ) CALL ADAMS_BASHFORTH3(
290           kkey = (itdkey-1)*Nr + k       I                                  bi, bj, 0,
291  #endif /* ALLOW_AUTODIFF_TAMC */       U                                  theta, gtNm,
292         I                                  tempStartAB, iterNb, myThid )
293  C--       Integrate continuity vertically for vertical velocity  C       compute S^n+1/2 (stored in gsNm) extrapolating S forward in time
294  c         CALL INTEGRATE_FOR_W(          IF ( AdamsBashforth_S ) CALL ADAMS_BASHFORTH3(
295  c    I                         bi, bj, k, uVel, vVel,       I                                  bi, bj, 0,
296  c    O                         wVel,       U                                  salt, gsNm,
297  c    I                         myThid )       I                                  saltStartAB, iterNb, myThid )
298    #endif /* ALLOW_ADAMSBASHFORTH_3 */
299  #ifdef    ALLOW_OBCS  
300  #ifdef    ALLOW_NONHYDROSTATIC  c       iMin = 1-OLx
301  C--       Apply OBC to W if in N-H mode  c       iMax = sNx+OLx
302  c         IF (useOBCS.AND.nonHydrostatic) THEN  c       jMin = 1-OLy
303  c           CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )  c       jMax = sNy+OLy
 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  
 #ifndef DISABLE_DEBUGMODE  
             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  
 #ifndef DISABLE_DEBUGMODE  
             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  
 #ifndef DISABLE_DEBUGMODE  
             IF ( debugLevel .GE. debLevB )  
      &       CALL DEBUG_CALL('CALC_IVDC',myThid)  
 #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 */  
   
 C--     end of diagnostic k loop (Nr:1)  
         ENDDO  
304    
305  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
306  cph avoids recomputation of integrate_for_w  cph avoids recomputation of integrate_for_w
307  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
308  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
309    
 #ifdef  ALLOW_OBCS  
 C--     Calculate future values on open boundaries  
         IF (useOBCS) THEN  
 #ifndef DISABLE_DEBUGMODE  
           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 */  
   
   
 #ifdef ALLOW_THERM_SEAICE  
        IF (useThermSeaIce) THEN  
 #ifndef DISABLE_DEBUGMODE  
         IF ( debugLevel .GE. debLevB )  
      &    CALL DEBUG_CALL('ICE_FORCING',myThid)  
 #endif  
 C--     Determines forcing terms based on external fields  
 C       including effects from ice  
         CALL ICE_FORCING(  
      I             bi, bj, iMin, iMax, jMin, jMax,  
      I             myThid )  
        ELSE  
 #endif /* ALLOW_THERM_SEAICE */  
   
 C--     Determines forcing terms based on external fields  
 C       relaxation terms, etc.  
 #ifndef DISABLE_DEBUGMODE  
         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_THERM_SEAICE  
 C--    end of if/else block useThermSeaIce --  
        ENDIF  
 #endif /* ALLOW_THERM_SEAICE */  
   
 #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  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
310  C--     Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h  C--     Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h
311  C--     MOST of THERMODYNAMICS will be disabled  C--     MOST of THERMODYNAMICS will be disabled
312  #ifndef SINGLE_LAYER_MODE  #ifndef SINGLE_LAYER_MODE
313    
 #ifdef  ALLOW_GMREDI  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 cph storing here is needed only for one GMREDI_OPTIONS:  
 cph define GM_BOLUS_ADVEC  
 cph but I've avoided the #ifdef for now, in case more things change  
 CADJ STORE sigmaX(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte  
 CADJ STORE sigmaY(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte  
 CADJ STORE sigmaR(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 C--     Calculate iso-neutral slopes for the GM/Redi parameterisation  
         IF (useGMRedi) THEN  
 #ifndef DISABLE_DEBUGMODE  
           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  
 #ifndef DISABLE_DEBUGMODE  
           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  
   
314  #ifdef ALLOW_AUTODIFF_TAMC  #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 */  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE KappaRT(:,:,:)     = comlev1_bibj, key=itdkey, byte=isbyte  
 CADJ STORE KappaRS(:,:,:)     = comlev1_bibj, key=itdkey, byte=isbyte  
315  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
316  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
317  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
318  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
319  #ifdef ALLOW_PASSIVE_TRACER  # ifdef ALLOW_DEPTH_CONTROL
320  CADJ STORE tr1  (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE gtnm1(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
321  #endif  CADJ STORE gsnm1(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
322  #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  
323  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
324    
 #ifdef ALLOW_AIM  
 C       AIM - atmospheric intermediate model, physics package code.  
         IF ( useAIM ) THEN  
 #ifndef DISABLE_DEBUGMODE  
           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 */  
   
325  #ifndef DISABLE_MULTIDIM_ADVECTION  #ifndef DISABLE_MULTIDIM_ADVECTION
326  C--     Some advection schemes are better calculated using a multi-dimensional  C--     Some advection schemes are better calculated using a multi-dimensional
327  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 566  C recomputation. It *is* differentiable, Line 334  C recomputation. It *is* differentiable,
334  C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to  C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to
335  C disable this section of code.  C disable this section of code.
336          IF (tempMultiDimAdvec) THEN          IF (tempMultiDimAdvec) THEN
337  #ifndef DISABLE_DEBUGMODE  #ifdef ALLOW_DEBUG
338            IF ( debugLevel .GE. debLevB )            IF ( debugLevel .GE. debLevB )
339       &     CALL DEBUG_CALL('GAD_ADVECTION',myThid)       &     CALL DEBUG_CALL('GAD_ADVECTION',myThid)
340  #endif  #endif
341            CALL GAD_ADVECTION(bi,bj,tempAdvScheme,GAD_TEMPERATURE,            CALL GAD_ADVECTION(
342       U                      theta,gT,       I             tempImplVertAdv, tempAdvScheme, tempVertAdvScheme,
343       I                      myTime,myIter,myThid)       I             GAD_TEMPERATURE,
344         I             uVel, vVel, wVel, theta,
345         O             gT,
346         I             bi,bj,myTime,myIter,myThid)
347          ENDIF          ENDIF
348          IF (saltMultiDimAdvec) THEN          IF (saltMultiDimAdvec) THEN
349  #ifndef DISABLE_DEBUGMODE  #ifdef ALLOW_DEBUG
350            IF ( debugLevel .GE. debLevB )            IF ( debugLevel .GE. debLevB )
351       &     CALL DEBUG_CALL('GAD_ADVECTION',myThid)       &     CALL DEBUG_CALL('GAD_ADVECTION',myThid)
352  #endif  #endif
353            CALL GAD_ADVECTION(bi,bj,saltAdvScheme,GAD_SALINITY,            CALL GAD_ADVECTION(
354       U                      salt,gS,       I             saltImplVertAdv, saltAdvScheme, saltVertAdvScheme,
355       I                      myTime,myIter,myThid)       I             GAD_SALINITY,
356         I             uVel, vVel, wVel, salt,
357         O             gS,
358         I             bi,bj,myTime,myIter,myThid)
359          ENDIF          ENDIF
360    
361  C Since passive tracers are configurable separately from T,S we  C Since passive tracers are configurable separately from T,S we
362  C call the multi-dimensional method for PTRACERS regardless  C call the multi-dimensional method for PTRACERS regardless
363  C of whether multiDimAdvection is set or not.  C of whether multiDimAdvection is set or not.
364  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
365          IF ( usePTRACERS ) THEN          IF ( usePTRACERS ) THEN
366  #ifndef DISABLE_DEBUGMODE  #ifdef ALLOW_DEBUG
367            IF ( debugLevel .GE. debLevB )            IF ( debugLevel .GE. debLevB )
368       &     CALL DEBUG_CALL('PTRACERS_ADVECTION',myThid)       &     CALL DEBUG_CALL('PTRACERS_ADVECTION',myThid)
369  #endif  #endif
# Line 597  C of whether multiDimAdvection is set or Line 372  C of whether multiDimAdvection is set or
372  #endif /* ALLOW_PTRACERS */  #endif /* ALLOW_PTRACERS */
373  #endif /* DISABLE_MULTIDIM_ADVECTION */  #endif /* DISABLE_MULTIDIM_ADVECTION */
374    
375  #ifndef DISABLE_DEBUGMODE  #ifdef ALLOW_DEBUG
376         IF ( debugLevel .GE. debLevB )          IF ( debugLevel .GE. debLevB )
377       &    CALL DEBUG_MSG('ENTERING DOWNWARD K LOOP',myThid)       &    CALL DEBUG_MSG('ENTERING DOWNWARD K LOOP',myThid)
378  #endif  #endif
379    
# Line 607  C--     Start of thermodynamics loop Line 382  C--     Start of thermodynamics loop
382  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
383  C? Patrick Is this formula correct?  C? Patrick Is this formula correct?
384  cph Yes, but I rewrote it.  cph Yes, but I rewrote it.
385  cph Also, the KappaR? need the index and subscript k!  cph Also, the kappaR? need the index and subscript k!
386           kkey = (itdkey-1)*Nr + k           kkey = (itdkey-1)*Nr + k
387  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
388    
# Line 624  C--       kDown  Cycles through 2,1 to p Line 399  C--       kDown  Cycles through 2,1 to p
399            jMin = 1-OLy            jMin = 1-OLy
400            jMax = sNy+OLy            jMax = sNy+OLy
401    
402  C--       Get temporary terms used by tendency routines            IF (k.EQ.Nr) THEN
403               DO j=1-Oly,sNy+Oly
404                DO i=1-Olx,sNx+Olx
405                 rTransKp1(i,j) = 0. _d 0
406                ENDDO
407               ENDDO
408              ELSE
409               DO j=1-Oly,sNy+Oly
410                DO i=1-Olx,sNx+Olx
411                 rTransKp1(i,j) = rTrans(i,j)
412                ENDDO
413               ENDDO
414              ENDIF
415    #ifdef ALLOW_AUTODIFF_TAMC
416    CADJ STORE rTransKp1(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
417    #endif
418    
419    C--       Get temporary terms used by tendency routines :
420    C-        Calculate horizontal "volume transport" through tracer cell face
421    C         anelastic: uTrans,vTrans are scaled by rhoFacC (~ mass transport)
422            CALL CALC_COMMON_FACTORS (            CALL CALC_COMMON_FACTORS (
423       I         bi,bj,iMin,iMax,jMin,jMax,k,       I         uVel, vVel,
424       O         xA,yA,uTrans,vTrans,rTrans,maskUp,       O         uFld, vFld, uTrans, vTrans, xA, yA,
425       I         myThid)       I         k,bi,bj, myThid )
426    
427    C-        Calculate vertical "volume transport" through tracer cell face
428              IF (k.EQ.1) THEN
429    C-         Surface interface :
430               DO j=1-Oly,sNy+Oly
431                DO i=1-Olx,sNx+Olx
432                 wFld(i,j)   = 0. _d 0
433                 maskUp(i,j) = 0. _d 0
434                 rTrans(i,j) = 0. _d 0
435                ENDDO
436               ENDDO
437              ELSE
438    C-         Interior interface :
439    C          anelastic: rTrans is scaled by rhoFacF (~ mass transport)
440               DO j=1-Oly,sNy+Oly
441                DO i=1-Olx,sNx+Olx
442                 wFld(i,j)   = wVel(i,j,k,bi,bj)
443                 maskUp(i,j) = maskC(i,j,k-1,bi,bj)*maskC(i,j,k,bi,bj)
444                 rTrans(i,j) = wFld(i,j)*rA(i,j,bi,bj)*maskUp(i,j)
445         &                              *deepFac2F(k)*rhoFacF(k)
446                ENDDO
447               ENDDO
448              ENDIF
449    
450  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
   
451  C--   Residual transp = Bolus transp + Eulerian transp  C--   Residual transp = Bolus transp + Eulerian transp
452            IF (useGMRedi) THEN            IF (useGMRedi) THEN
453              CALL GMREDI_CALC_UVFLOW(              CALL GMREDI_CALC_UVFLOW(
454       &            uTrans, vTrans, bi, bj, k, myThid)       U                  uFld, vFld, uTrans, vTrans,
455              IF (K.GE.2) CALL GMREDI_CALC_WFLOW(       I                  k, bi, bj, myThid )
456       &                    rTrans, bi, bj, k, myThid)              IF (K.GE.2) THEN
457                  CALL GMREDI_CALC_WFLOW(
458         U                  wFld, rTrans,
459         I                  k, bi, bj, myThid )
460                ENDIF
461            ENDIF            ENDIF
462    # ifdef ALLOW_AUTODIFF_TAMC
463  #ifdef ALLOW_AUTODIFF_TAMC  CADJ STORE rTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
464  #ifdef GM_BOLUS_ADVEC  CADJ STORE wfld(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
465    # ifdef GM_BOLUS_ADVEC
466    CADJ STORE ufld(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
467    CADJ STORE vfld(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
468  CADJ STORE uTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE uTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
469  CADJ STORE vTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE vTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
470  CADJ STORE rTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte  # endif
471  #endif  # endif /* ALLOW_AUTODIFF_TAMC */
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
472  #endif /* ALLOW_GMREDI */  #endif /* ALLOW_GMREDI */
473    
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE KappaRT(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte  
 CADJ STORE KappaRS(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
474  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL
475  C--      Calculate the total vertical diffusivity  C--      Calculate the total vertical diffusivity
476           CALL CALC_DIFFUSIVITY(            IF ( .NOT.implicitDiffusion ) THEN
477       I        bi,bj,iMin,iMax,jMin,jMax,k,              CALL CALC_DIFFUSIVITY(
478       I        maskUp,       I          bi,bj,iMin,iMax,jMin,jMax,k,
479       O        KappaRT,KappaRS,       I          maskUp,
480       I        myThid)       O          kappaRT,kappaRS,
481         I          myThid)
482              ENDIF
483    # ifdef ALLOW_AUTODIFF_TAMC
484    CADJ STORE kappaRT(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
485    CADJ STORE kappaRS(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
486    # endif /* ALLOW_AUTODIFF_TAMC */
487  #endif  #endif
488    
489            iMin = 1-OLx+2            iMin = 1-OLx+2
# Line 670  C--      Calculate the total vertical di Line 492  C--      Calculate the total vertical di
492            jMax = sNy+OLy-1            jMax = sNy+OLy-1
493    
494  C--      Calculate active tracer tendencies (gT,gS,...)  C--      Calculate active tracer tendencies (gT,gS,...)
495  C        and step forward storing result in gTnm1, gSnm1, etc.  C        and step forward storing result in gT, gS, etc.
496    C--
497    # ifdef ALLOW_AUTODIFF_TAMC
498    #  if ((defined NONLIN_FRSURF) || (defined ALLOW_DEPTH_CONTROL)) && (defined ALLOW_GMREDI)
499    CADJ STORE kux(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
500    CADJ STORE kvy(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
501    #   ifdef GM_EXTRA_DIAGONAL
502    CADJ STORE kuz(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
503    CADJ STORE kvz(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
504    #   endif
505    #  endif
506    # endif /* ALLOW_AUTODIFF_TAMC */
507    C
508    #ifdef ALLOW_AUTODIFF_TAMC
509    # if (defined NONLIN_FRSURF) || (defined ALLOW_DEPTH_CONTROL)
510    cph-test
511    CADJ STORE uFld(:,:), vFld(:,:), wFld(:,:)
512    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
513    CADJ STORE uTrans(:,:), vTrans(:,:)
514    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
515    CADJ STORE xA(:,:), yA(:,:)
516    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
517    # endif
518    #endif /* ALLOW_AUTODIFF_TAMC */
519    C
520           IF ( tempStepping ) THEN           IF ( tempStepping ) THEN
521             CALL CALC_GT(  #ifdef ALLOW_AUTODIFF_TAMC
522    CADJ STORE gTnm1(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
523    # if (defined NONLIN_FRSURF) || (defined ALLOW_DEPTH_CONTROL)
524    CADJ STORE fvert(:,:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
525    # endif
526    #endif /* ALLOW_AUTODIFF_TAMC */
527              CALL CALC_GT(
528       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
529       I         xA,yA,uTrans,vTrans,rTrans,maskUp,       I         xA, yA, maskUp, uFld, vFld, wFld,
530       I         KappaRT,       I         uTrans, vTrans, rTrans, rTransKp1,
531         I         kappaRT,
532       U         fVerT,       U         fVerT,
533       I         myTime,myIter,myThid)       I         myTime,myIter,myThid)
534    #ifdef ALLOW_ADAMSBASHFORTH_3
535              IF ( AdamsBashforth_T ) THEN
536             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
537       I         bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,       I         bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,
538       I         theta, gT,       I         gtNm(1-Olx,1-Oly,1,1,1,m2),
539         U         gT,
540       I         myIter, myThid)       I         myIter, myThid)
541           ENDIF            ELSE
542    #endif
543  #ifdef ALLOW_THERM_SEAICE             CALL TIMESTEP_TRACER(
544           IF (useThermSeaIce .AND. k.EQ.1) THEN       I         bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,
545             CALL ICE_FREEZE( bi,bj, iMin,iMax,jMin,jMax, myThid )       I         theta,
546           ENDIF       U         gT,
547         I         myIter, myThid)
548    #ifdef ALLOW_ADAMSBASHFORTH_3
549              ENDIF
550  #endif  #endif
551             ENDIF
552    
553           IF ( saltStepping ) THEN           IF ( saltStepping ) THEN
554             CALL CALC_GS(  #ifdef ALLOW_AUTODIFF_TAMC
555    CADJ STORE gSnm1(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
556    # if (defined NONLIN_FRSURF) || (defined ALLOW_DEPTH_CONTROL)
557    CADJ STORE fvers(:,:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
558    # endif
559    #endif /* ALLOW_AUTODIFF_TAMC */
560              CALL CALC_GS(
561       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
562       I         xA,yA,uTrans,vTrans,rTrans,maskUp,       I         xA, yA, maskUp, uFld, vFld, wFld,
563       I         KappaRS,       I         uTrans, vTrans, rTrans, rTransKp1,
564         I         kappaRS,
565       U         fVerS,       U         fVerS,
566       I         myTime,myIter,myThid)       I         myTime,myIter,myThid)
567    #ifdef ALLOW_ADAMSBASHFORTH_3
568              IF ( AdamsBashforth_S ) THEN
569             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
570       I         bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,       I         bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,
571       I         salt, gS,       I         gsNm(1-Olx,1-Oly,1,1,1,m2),
572         U         gS,
573       I         myIter, myThid)       I         myIter, myThid)
574           ENDIF            ELSE
575  #ifdef ALLOW_PASSIVE_TRACER  #endif
          IF ( tr1Stepping ) THEN  
            CALL CALC_GTR1(  
      I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,  
      I         xA,yA,uTrans,vTrans,rTrans,maskUp,  
      I         KappaRT,  
      U         fVerTr1,  
      I         myTime,myIter,myThid)  
576             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
577       I         bi,bj,iMin,iMax,jMin,jMax,k,tracerAdvScheme,       I         bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,
578       I         Tr1, gTr1,       I         salt,
579       I         myIter,myThid)       U         gS,
580           ENDIF       I         myIter, myThid)
581    #ifdef ALLOW_ADAMSBASHFORTH_3
582              ENDIF
583  #endif  #endif
584             ENDIF
585    
586  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
587           IF ( usePTRACERS ) THEN           IF ( usePTRACERS ) THEN
588               IF ( .NOT.implicitDiffusion ) THEN
589                 CALL PTRACERS_CALC_DIFF(
590         I            bi,bj,iMin,iMax,jMin,jMax,k,
591         I            maskUp,
592         O            kappaRTr,
593         I            myThid)
594               ENDIF
595    # ifdef ALLOW_AUTODIFF_TAMC
596    CADJ STORE kappaRTr(:,:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
597    # endif /* ALLOW_AUTODIFF_TAMC */
598             CALL PTRACERS_INTEGRATE(             CALL PTRACERS_INTEGRATE(
599       I         bi,bj,k,       I          bi,bj,k,
600       I         xA,yA,uTrans,vTrans,rTrans,maskUp,       I          xA, yA, maskUp, uFld, vFld, wFld,
601       X         KappaRS,       I          uTrans, vTrans, rTrans, rTransKp1,
602       I         myIter,myTime,myThid)       I          kappaRTr,
603         U          fVerP,
604         I          myTime,myIter,myThid)
605           ENDIF           ENDIF
606  #endif /* ALLOW_PTRACERS */  #endif /* ALLOW_PTRACERS */
607    
# Line 734  C--      Apply open boundary conditions Line 613  C--      Apply open boundary conditions
613  #endif   /* ALLOW_OBCS */  #endif   /* ALLOW_OBCS */
614    
615  C--      Freeze water  C--      Freeze water
616           IF ( allowFreezing .AND. .NOT. useSEAICE  C  this bit of code is left here for backward compatibility.
617       &       .AND. .NOT.(useThermSeaIce.AND.k.EQ.1) ) THEN  C  freezing at surface level has been moved to FORWARD_STEP
618             IF ( useOldFreezing .AND. .NOT. useSEAICE
619         &       .AND. .NOT.(useThSIce.AND.k.EQ.1) ) THEN
620  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
621  CADJ STORE gT(:,:,k,bi,bj) = comlev1_bibj_k  CADJ STORE gT(:,:,k,bi,bj) = comlev1_bibj_k
622  CADJ &   , key = kkey, byte = isbyte  CADJ &   , key = kkey, byte = isbyte
# Line 746  CADJ &   , key = kkey, byte = isbyte Line 627  CADJ &   , key = kkey, byte = isbyte
627  C--     end of thermodynamic k loop (Nr:1)  C--     end of thermodynamic k loop (Nr:1)
628          ENDDO          ENDDO
629    
630  cswdice -- add ---  C       All explicit advection/diffusion/sources should now be
631  #ifdef ALLOW_THERM_SEAICE  C       done. The updated tracer field is in gPtr. Accumalate
632  c timeaveraging for ice model values  C       explicit tendency and also reset gPtr to initial tracer
633             CALL ICE_AVE(bi,bj,iMin,iMax,jMin,jMax,myThid )  C       field for implicit matrix calculation
634  #endif  
635  cswdice --- end add ---  #ifdef ALLOW_MATRIX
636            IF (useMATRIX)
637         &    CALL MATRIX_STORE_TENDENCY_EXP(bi,bj, myTime,myIter,myThid)
638    #endif
639    
640  C--     Implicit diffusion          iMin = 1
641          IF (implicitDiffusion) THEN          iMax = sNx
642            jMin = 1
643           IF (tempStepping) THEN          jMax = sNy
644    
645    C--     Implicit vertical advection & diffusion
646            IF ( tempStepping .AND. implicitDiffusion ) THEN
647              CALL CALC_3D_DIFFUSIVITY(
648         I         bi,bj,iMin,iMax,jMin,jMax,
649         I         GAD_TEMPERATURE, useGMredi, useKPP,
650         O         kappaRk,
651         I         myThid)
652            ENDIF
653    #ifdef INCLUDE_IMPLVERTADV_CODE
654            IF ( tempImplVertAdv ) THEN
655              CALL GAD_IMPLICIT_R(
656         I         tempImplVertAdv, tempAdvScheme, GAD_TEMPERATURE,
657         I         kappaRk, wVel, theta,
658         U         gT,
659         I         bi, bj, myTime, myIter, myThid )
660            ELSEIF ( tempStepping .AND. implicitDiffusion ) THEN
661    #else /* INCLUDE_IMPLVERTADV_CODE */
662            IF     ( tempStepping .AND. implicitDiffusion ) THEN
663    #endif /* INCLUDE_IMPLVERTADV_CODE */
664  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
665    CADJ STORE kappaRk(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
666  CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte  CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
667  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
668              CALL IMPLDIFF(            CALL IMPLDIFF(
669       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
670       I         deltaTtracer, KappaRT, recip_HFacC,       I         GAD_TEMPERATURE, kappaRk, recip_hFacC,
671       U         gT,       U         gT,
672       I         myThid )       I         myThid )
673            ENDIF
674    
675    #ifdef ALLOW_TIMEAVE
676            useVariableK = useKPP .OR. usePP81 .OR. useMY82 .OR. useGGL90
677         &       .OR. useGMredi .OR. ivdc_kappa.NE.0.
678            IF (taveFreq.GT.0. .AND. useVariableK ) THEN
679             IF (implicitDiffusion) THEN
680               CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, gT, kappaRk,
681         I                        Nr, 3, deltaTclock, bi, bj, myThid)
682    c        ELSE
683    c          CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, theta, kappaRT,
684    c    I                        Nr, 3, deltaTclock, bi, bj, myThid)
685           ENDIF           ENDIF
686            ENDIF
687    #endif /* ALLOW_TIMEAVE */
688    
689            IF ( saltStepping .AND. implicitDiffusion ) THEN
690              CALL CALC_3D_DIFFUSIVITY(
691         I         bi,bj,iMin,iMax,jMin,jMax,
692         I         GAD_SALINITY, useGMredi, useKPP,
693         O         kappaRk,
694         I         myThid)
695            ENDIF
696    
697           IF (saltStepping) THEN  #ifdef INCLUDE_IMPLVERTADV_CODE
698            IF ( saltImplVertAdv ) THEN
699              CALL GAD_IMPLICIT_R(
700         I         saltImplVertAdv, saltAdvScheme, GAD_SALINITY,
701         I         kappaRk, wVel, salt,
702         U         gS,
703         I         bi, bj, myTime, myIter, myThid )
704            ELSEIF ( saltStepping .AND. implicitDiffusion ) THEN
705    #else /* INCLUDE_IMPLVERTADV_CODE */
706            IF     ( saltStepping .AND. implicitDiffusion ) THEN
707    #endif /* INCLUDE_IMPLVERTADV_CODE */
708  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
709    CADJ STORE kappaRk(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
710  CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte  CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
711  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
712              CALL IMPLDIFF(            CALL IMPLDIFF(
713       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
714       I         deltaTtracer, KappaRS, recip_HFacC,       I         GAD_SALINITY, kappaRk, recip_hFacC,
715       U         gS,       U         gS,
716       I         myThid )       I         myThid )
717           ENDIF          ENDIF
   
 #ifdef ALLOW_PASSIVE_TRACER  
          IF (tr1Stepping) 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  
718    
719  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
720  C Vertical diffusion (implicit) for passive tracers          IF     ( usePTRACERS ) THEN
721           IF ( usePTRACERS ) THEN  C--     Vertical advection/diffusion (implicit) for passive tracers
722             CALL PTRACERS_IMPLDIFF( bi,bj,KappaRS,myThid )             CALL PTRACERS_IMPLICIT(
723           ENDIF       U                             kappaRk,
724         I                             bi, bj, myTime, myIter, myThid )
725            ENDIF
726  #endif /* ALLOW_PTRACERS */  #endif /* ALLOW_PTRACERS */
727    
728  #ifdef   ALLOW_OBCS  #ifdef   ALLOW_OBCS
729  C--      Apply open boundary conditions  C--      Apply open boundary conditions
730           IF (useOBCS) THEN          IF ( ( implicitDiffusion
731         &    .OR. tempImplVertAdv
732         &    .OR. saltImplVertAdv
733         &       ) .AND. useOBCS     ) THEN
734             DO K=1,Nr             DO K=1,Nr
735               CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )               CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )
736             ENDDO             ENDDO
          END IF  
 #endif   /* ALLOW_OBCS */  
   
 C--     End If implicitDiffusion  
737          ENDIF          ENDIF
738    #endif   /* ALLOW_OBCS */
 #ifdef ALLOW_TIMEAVE  
         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 /* ALLOW_TIMEAVE */  
739    
740  #endif /* SINGLE_LAYER_MODE */  #endif /* SINGLE_LAYER_MODE */
741    
# Line 836  C--   end bi,bj loops. Line 743  C--   end bi,bj loops.
743         ENDDO         ENDDO
744        ENDDO        ENDDO
745    
746  #ifndef DISABLE_DEBUGMODE  #ifdef ALLOW_DEBUG
747        If (debugMode) THEN        If (debugMode) THEN
748         CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (THERMODYNAMICS)',myThid)
749         CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (THERMODYNAMICS)',myThid)
750         CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (THERMODYNAMICS)',myThid)
751         CALL DEBUG_STATS_RL(Nr,theta,'Theta (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,theta,'Theta (THERMODYNAMICS)',myThid)
752         CALL DEBUG_STATS_RL(Nr,salt,'Salt (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,salt,'Salt (THERMODYNAMICS)',myThid)
753         CALL DEBUG_STATS_RL(Nr,Gt,'Gt (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,gT,'Gt (THERMODYNAMICS)',myThid)
754         CALL DEBUG_STATS_RL(Nr,Gs,'Gs (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,gS,'Gs (THERMODYNAMICS)',myThid)
755         CALL DEBUG_STATS_RL(Nr,GtNm1,'GtNm1 (THERMODYNAMICS)',myThid)  #ifndef ALLOW_ADAMSBASHFORTH_3
756         CALL DEBUG_STATS_RL(Nr,GsNm1,'GsNm1 (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,gtNm1,'GtNm1 (THERMODYNAMICS)',myThid)
757           CALL DEBUG_STATS_RL(Nr,gsNm1,'GsNm1 (THERMODYNAMICS)',myThid)
758    #endif
759  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
760         IF ( usePTRACERS ) THEN         IF ( usePTRACERS ) THEN
761           CALL PTRACERS_DEBUG(myThid)           CALL PTRACERS_DEBUG(myThid)
762         ENDIF         ENDIF
763  #endif /* ALLOW_PTRACERS */  #endif /* ALLOW_PTRACERS */
764        ENDIF        ENDIF
765  #endif  #endif /* ALLOW_DEBUG */
766    
767  #ifndef DISABLE_DEBUGMODE  #ifdef ALLOW_DEBUG
768           IF ( debugLevel .GE. debLevB )           IF ( debugLevel .GE. debLevB )
769       &    CALL DEBUG_LEAVE('FORWARD_STEP',myThid)       &    CALL DEBUG_LEAVE('THERMODYNAMICS',myThid)
770  #endif  #endif
771    
772    #endif /* ALLOW_GENERIC_ADVDIFF */
773    
774        RETURN        RETURN
775        END        END
Legend:
Removed from v.1.50  
changed lines
  Added in v.1.109
  ViewVC Help
Powered by ViewVC 1.1.22