/[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.14 by jmc, Fri Nov 16 03:25:41 2001 UTC revision 1.72 by jmc, Tue Jul 6 00:58:40 2004 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    #ifdef ALLOW_PTRACERS
7    # include "PTRACERS_OPTIONS.h"
8    #endif
9    
10    #ifdef ALLOW_AUTODIFF_TAMC
11    # ifdef ALLOW_GMREDI
12    #  include "GMREDI_OPTIONS.h"
13    # endif
14    # ifdef ALLOW_KPP
15    #  include "KPP_OPTIONS.h"
16    # endif
17    #endif /* ALLOW_AUTODIFF_TAMC */
18    
19  CBOP  CBOP
20  C     !ROUTINE: THERMODYNAMICS  C     !ROUTINE: THERMODYNAMICS
# Line 71  C     == Global variables === Line 84  C     == Global variables ===
84  #ifdef ALLOW_PASSIVE_TRACER  #ifdef ALLOW_PASSIVE_TRACER
85  #include "TR1.h"  #include "TR1.h"
86  #endif  #endif
87    #ifdef ALLOW_PTRACERS
88    #include "PTRACERS.h"
89    #endif
90    #ifdef ALLOW_TIMEAVE
91    #include "TIMEAVE_STATV.h"
92    #endif
93    
94  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
95  # include "tamc.h"  # include "tamc.h"
96  # include "tamc_keys.h"  # include "tamc_keys.h"
97  # include "FFIELDS.h"  # include "FFIELDS.h"
98    # include "EOS.h"
99  # ifdef ALLOW_KPP  # ifdef ALLOW_KPP
100  #  include "KPP.h"  #  include "KPP.h"
101  # endif  # endif
102  # ifdef ALLOW_GMREDI  # ifdef ALLOW_GMREDI
103  #  include "GMREDI.h"  #  include "GMREDI.h"
104  # endif  # endif
105    # ifdef ALLOW_EBM
106    #  include "EBM.h"
107    # endif
108  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
 #ifdef ALLOW_TIMEAVE  
 #include "TIMEAVE_STATV.h"  
 #endif  
109    
110  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
111  C     == Routine arguments ==  C     == Routine arguments ==
# Line 103  C                              transport Line 124  C                              transport
124  C                              o uTrans: Zonal transport  C                              o uTrans: Zonal transport
125  C                              o vTrans: Meridional transport  C                              o vTrans: Meridional transport
126  C                              o rTrans: Vertical transport  C                              o rTrans: Vertical transport
127    C     rTransKp1                o vertical volume transp. at interface k+1
128  C     maskUp                   o maskUp: land/water mask for W points  C     maskUp                   o maskUp: land/water mask for W points
129  C     fVer[STUV]               o fVer: Vertical flux term - note fVer  C     fVer[STUV]               o fVer: Vertical flux term - note fVer
130  C                                      is "pipelined" in the vertical  C                                      is "pipelined" in the vertical
131  C                                      so we need an fVer for each  C                                      so we need an fVer for each
132  C                                      variable.  C                                      variable.
 C     rhoK, rhoKM1   - Density at current level, and level above  
 C     phiHyd         - Hydrostatic part of the potential phiHydi.  
 C                      In z coords phiHydiHyd is the hydrostatic  
 C                      Potential (=pressure/rho0) anomaly  
 C                      In p coords phiHydiHyd is the geopotential  
 C                      surface height anomaly.  
 C     phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean)  
 C     phiSurfY             or geopotentiel (atmos) in X and Y direction  
133  C     KappaRT,       - Total diffusion in vertical for T and S.  C     KappaRT,       - Total diffusion in vertical for T and S.
134  C     KappaRS          (background + spatially varying, isopycnal term).  C     KappaRS          (background + spatially varying, isopycnal term).
135    C     useVariableK   = T when vertical diffusion is not constant
136  C     iMin, iMax     - Ranges and sub-block indices on which calculations  C     iMin, iMax     - Ranges and sub-block indices on which calculations
137  C     jMin, jMax       are applied.  C     jMin, jMax       are applied.
138  C     bi, bj  C     bi, bj
# Line 129  C                      index into fVerTe Line 144  C                      index into fVerTe
144        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
145        _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
146        _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
147          _RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
148        _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
149        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
150        _RL fVerS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
151    #ifdef ALLOW_PASSIVE_TRACER
152        _RL fVerTr1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerTr1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
153        _RL phiHyd  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  #endif
154        _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  #ifdef ALLOW_PTRACERS
155        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fVerP   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2,PTRACERS_num)
156        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  #endif
       _RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
157        _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
158        _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
159        _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
160        _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
161        _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
162  C     This is currently used by IVDC and Diagnostics        _RL kp1Msk
163        _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        LOGICAL useVariableK
164        INTEGER iMin, iMax        INTEGER iMin, iMax
165        INTEGER jMin, jMax        INTEGER jMin, jMax
166        INTEGER bi, bj        INTEGER bi, bj
167        INTEGER i, j        INTEGER i, j
168        INTEGER k, km1, kup, kDown        INTEGER k, km1, kup, kDown
169          INTEGER iTracer, ip
170    
171  CEOP  CEOP
172    
173    #ifdef ALLOW_DEBUG
174             IF ( debugLevel .GE. debLevB )
175         &    CALL DEBUG_ENTER('THERMODYNAMICS',myThid)
176    #endif
177    
178  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
179  C--   dummy statement to end declaration part  C--   dummy statement to end declaration part
180        ikey = 1        ikey = 1
181          itdkey = 1
182  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
183    
 C--   Set up work arrays with valid (i.e. not NaN) values  
 C     These inital values do not alter the numerical results. They  
 C     just ensure that all memory references are to valid floating  
 C     point numbers. This prevents spurious hardware signals due to  
 C     uninitialised but inert locations.  
       DO j=1-OLy,sNy+OLy  
        DO i=1-OLx,sNx+OLx  
         xA(i,j)      = 0. _d 0  
         yA(i,j)      = 0. _d 0  
         uTrans(i,j)  = 0. _d 0  
         vTrans(i,j)  = 0. _d 0  
         DO k=1,Nr  
          phiHyd(i,j,k)  = 0. _d 0  
          sigmaX(i,j,k) = 0. _d 0  
          sigmaY(i,j,k) = 0. _d 0  
          sigmaR(i,j,k) = 0. _d 0  
         ENDDO  
         rhoKM1 (i,j) = 0. _d 0  
         rhok   (i,j) = 0. _d 0  
         phiSurfX(i,j) = 0. _d 0  
         phiSurfY(i,j) = 0. _d 0  
        ENDDO  
       ENDDO  
   
   
184  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
185  C--   HPF directive to help TAMC  C--   HPF directive to help TAMC
186  CHPF$ INDEPENDENT  CHPF$ INDEPENDENT
# Line 193  CHPF$ INDEPENDENT Line 191  CHPF$ INDEPENDENT
191  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
192  C--    HPF directive to help TAMC  C--    HPF directive to help TAMC
193  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS
194  CHPF$&                  ,phiHyd,utrans,vtrans,xA,yA  CHPF$&                  ,utrans,vtrans,xA,yA
195  CHPF$&                  ,KappaRT,KappaRS  CHPF$&                  ,KappaRT,KappaRS
196  CHPF$&                  )  CHPF$&                  )
197  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
# Line 208  CHPF$&                  ) Line 206  CHPF$&                  )
206            act3 = myThid - 1            act3 = myThid - 1
207            max3 = nTx*nTy            max3 = nTx*nTy
208            act4 = ikey_dynamics - 1            act4 = ikey_dynamics - 1
209            ikey = (act1 + 1) + act2*max1            itdkey = (act1 + 1) + act2*max1
210       &                      + act3*max1*max2       &                      + act3*max1*max2
211       &                      + act4*max1*max2*max3       &                      + act4*max1*max2*max3
212  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
213    
214  C--     Set up work arrays that need valid initial values  C--   Set up work arrays with valid (i.e. not NaN) values
215    C     These inital values do not alter the numerical results. They
216    C     just ensure that all memory references are to valid floating
217    C     point numbers. This prevents spurious hardware signals due to
218    C     uninitialised but inert locations.
219    
220          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
221           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
222              xA(i,j)        = 0. _d 0
223              yA(i,j)        = 0. _d 0
224              uTrans(i,j)    = 0. _d 0
225              vTrans(i,j)    = 0. _d 0
226            rTrans (i,j)   = 0. _d 0            rTrans (i,j)   = 0. _d 0
227              rTransKp1(i,j) = 0. _d 0
228            fVerT  (i,j,1) = 0. _d 0            fVerT  (i,j,1) = 0. _d 0
229            fVerT  (i,j,2) = 0. _d 0            fVerT  (i,j,2) = 0. _d 0
230            fVerS  (i,j,1) = 0. _d 0            fVerS  (i,j,1) = 0. _d 0
231            fVerS  (i,j,2) = 0. _d 0            fVerS  (i,j,2) = 0. _d 0
232    #ifdef ALLOW_PASSIVE_TRACER
233            fVerTr1(i,j,1) = 0. _d 0            fVerTr1(i,j,1) = 0. _d 0
234            fVerTr1(i,j,2) = 0. _d 0            fVerTr1(i,j,2) = 0. _d 0
235    #endif
236    #ifdef ALLOW_PTRACERS
237              DO ip=1,PTRACERS_num
238               fVerP  (i,j,1,ip) = 0. _d 0
239               fVerP  (i,j,2,ip) = 0. _d 0
240              ENDDO
241    #endif
242           ENDDO           ENDDO
243          ENDDO          ENDDO
244    
# Line 230  C--     Set up work arrays that need val Line 246  C--     Set up work arrays that need val
246           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
247            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
248  C This is currently also used by IVDC and Diagnostics  C This is currently also used by IVDC and Diagnostics
249             ConvectCount(i,j,k) = 0.             KappaRT(i,j,k)    = 0. _d 0
250             KappaRT(i,j,k) = 0. _d 0             KappaRS(i,j,k)    = 0. _d 0
251             KappaRS(i,j,k) = 0. _d 0  C- tracer tendency needs to be set to zero (moved here from gad_calc_rhs):
252  #ifdef ALLOW_AUTODIFF_TAMC             gT(i,j,k,bi,bj)   = 0. _d 0
253             gT(i,j,k,bi,bj) = 0. _d 0             gS(i,j,k,bi,bj)   = 0. _d 0
254             gS(i,j,k,bi,bj) = 0. _d 0  # ifdef ALLOW_PASSIVE_TRACER
255  #ifdef ALLOW_PASSIVE_TRACER  ceh3 needs an IF ( use PASSIVE_TRACER) THEN
256             gTr1(i,j,k,bi,bj) = 0. _d 0             gTr1(i,j,k,bi,bj) = 0. _d 0
257  #endif  # endif
258  #endif  # ifdef ALLOW_PTRACERS
259    ceh3 this should have an   IF ( usePTRACERS ) THEN
260               DO iTracer=1,PTRACERS_numInUse
261                gPTr(i,j,k,bi,bj,itracer) = 0. _d 0
262               ENDDO
263    # endif
264            ENDDO            ENDDO
265           ENDDO           ENDDO
266          ENDDO          ENDDO
267    
268          iMin = 1-OLx+1  c       iMin = 1-OLx
269          iMax = sNx+OLx  c       iMax = sNx+OLx
270          jMin = 1-OLy+1  c       jMin = 1-OLy
271          jMax = sNy+OLy  c       jMax = sNy+OLy
   
272    
273  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
274  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
275  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
276    CADJ STORE totphihyd
277    CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
278  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
279  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
280  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
281  #endif  #endif
282  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
283    
 C--     Start of diagnostic loop  
         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 = (ikey-1)*Nr + k  
 CADJ STORE rhokm1(:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte  
 CADJ STORE rhok  (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 C--       Integrate continuity vertically for vertical velocity  
           CALL INTEGRATE_FOR_W(  
      I                         bi, bj, k, uVel, vVel,  
      O                         wVel,  
      I                         myThid )  
   
 #ifdef    ALLOW_OBCS  
 #ifdef    ALLOW_NONHYDROSTATIC  
 C--       Apply OBC to W if in N-H mode  
           IF (useOBCS.AND.nonHydrostatic) THEN  
             CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )  
           ENDIF  
 #endif    /* ALLOW_NONHYDROSTATIC */  
 #endif    /* ALLOW_OBCS */  
   
 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_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, eosType,  
      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, eosType,  
      I        theta, salt,  
      O        rhoKm1,  
      I        myThid )  
             ENDIF  
             CALL GRAD_SIGMA(  
      I             bi, bj, iMin, iMax, jMin, jMax, k,  
      I             rhoK, rhoKm1, rhoK,  
      O             sigmaX, sigmaY, sigmaR,  
      I             myThid )  
           ENDIF  
   
 C--       Implicit Vertical Diffusion for Convection  
 c ==> should use sigmaR !!!  
           IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN  
             CALL CALC_IVDC(  
      I        bi, bj, iMin, iMax, jMin, jMax, k,  
      I        rhoKm1, rhoK,  
      U        ConvectCount, KappaRT, KappaRS,  
      I        myTime, myIter, myThid)  
           ENDIF  
   
 C--     end of diagnostic k loop (Nr:1)  
         ENDDO  
   
284  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
285  cph avoids recomputation of integrate_for_w  cph avoids recomputation of integrate_for_w
286  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
287  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
288    
289  #ifdef  ALLOW_OBCS  C--     Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h
290  C--     Calculate future values on open boundaries  C--     MOST of THERMODYNAMICS will be disabled
291          IF (useOBCS) THEN  #ifndef SINGLE_LAYER_MODE
           CALL OBCS_CALC( bi, bj, myTime+deltaT,  
      I            uVel, vVel, wVel, theta, salt,  
      I            myThid )  
         ENDIF  
 #endif  /* ALLOW_OBCS */  
   
 C--     Determines forcing terms based on external fields  
 C       relaxation terms, etc.  
         CALL EXTERNAL_FORCING_SURF(  
      I             bi, bj, iMin, iMax, jMin, jMax,  
      I             myThid )  
 #ifdef ALLOW_AUTODIFF_TAMC  
 cph needed for KPP  
 CADJ STORE surfacetendencyU(:,:,bi,bj)  
 CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE surfacetendencyV(:,:,bi,bj)  
 CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE surfacetendencyS(:,:,bi,bj)  
 CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE surfacetendencyT(:,:,bi,bj)  
 CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 #ifdef  ALLOW_GMREDI  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE sigmaX(:,:,:) = comlev1, key=ikey, byte=isbyte  
 CADJ STORE sigmaY(:,:,:) = comlev1, key=ikey, byte=isbyte  
 CADJ STORE sigmaR(:,:,:) = comlev1, key=ikey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
 C--     Calculate iso-neutral slopes for the GM/Redi parameterisation  
         IF (useGMRedi) THEN  
           DO k=1,Nr  
             CALL GMREDI_CALC_TENSOR(  
      I             bi, bj, iMin, iMax, jMin, jMax, k,  
      I             sigmaX, sigmaY, sigmaR,  
      I             myThid )  
           ENDDO  
 #ifdef ALLOW_AUTODIFF_TAMC  
         ELSE  
           DO k=1, Nr  
             CALL GMREDI_CALC_TENSOR_DUMMY(  
      I             bi, bj, iMin, iMax, jMin, jMax, k,  
      I             sigmaX, sigmaY, sigmaR,  
      I             myThid )  
           ENDDO  
 #endif /* ALLOW_AUTODIFF_TAMC */  
         ENDIF  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE Kwx(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE Kwy(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE Kwz(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 #endif  /* ALLOW_GMREDI */  
   
 #ifdef  ALLOW_KPP  
 C--     Compute KPP mixing coefficients  
         IF (useKPP) THEN  
           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  
292    
293  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
294  CADJ STORE KPPghat   (:,:,:,bi,bj)  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
295  CADJ &   , KPPdiffKzT(:,:,:,bi,bj)  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
296  CADJ &   , KPPdiffKzS(:,:,:,bi,bj)  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
297  CADJ &   , KPPfrac   (:,:  ,bi,bj)  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
 CADJ &                 = comlev1_bibj, key=ikey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 #endif  /* ALLOW_KPP */  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE KappaRT(:,:,:)     = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE KappaRS(:,:,:)     = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  
298  #ifdef ALLOW_PASSIVE_TRACER  #ifdef ALLOW_PASSIVE_TRACER
299  CADJ STORE tr1  (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE tr1  (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
300    #endif
301    #ifdef ALLOW_PTRACERS
302    cph-- moved to forward_step to avoid key computation
303    cphCADJ STORE ptracer(:,:,:,bi,bj,itracer) = comlev1_bibj,
304    cphCADJ &                              key=itdkey, byte=isbyte
305  #endif  #endif
306  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
307    
 #ifdef ALLOW_AIM  
 C       AIM - atmospheric intermediate model, physics package code.  
 C note(jmc) : phiHyd=0 at this point but is not really used in Molteni Physics  
         IF ( useAIM ) THEN  
          CALL TIMER_START('AIM_DO_ATMOS_PHYS      [DYNAMICS]', myThid)  
          CALL AIM_DO_ATMOS_PHYSICS( phiHyd, bi, bj, myTime, myThid )  
          CALL TIMER_STOP ('AIM_DO_ATMOS_PHYS      [DYNAMICS]', myThid)  
         ENDIF  
 #endif /* ALLOW_AIM */  
   
 #ifdef ALLOW_TIMEAVE  
         IF (taveFreq.GT.0. .AND. ivdc_kappa.NE.0.) THEN  
           CALL TIMEAVE_CUMULATE(ConvectCountTave, ConvectCount, Nr,  
      I                               deltaTclock, bi, bj, myThid)  
         ENDIF  
 #endif /* ALLOW_TIMEAVE */  
   
308  #ifndef DISABLE_MULTIDIM_ADVECTION  #ifndef DISABLE_MULTIDIM_ADVECTION
309  C--     Some advection schemes are better calculated using a multi-dimensional  C--     Some advection schemes are better calculated using a multi-dimensional
310  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 458  C to be able to exclude this scheme to a Line 316  C to be able to exclude this scheme to a
316  C recomputation. It *is* differentiable, if you need it.  C recomputation. It *is* differentiable, if you need it.
317  C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to  C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to
318  C disable this section of code.  C disable this section of code.
319          IF (multiDimAdvection) THEN          IF (tempMultiDimAdvec) THEN
320           IF (tempStepping .AND.  #ifdef ALLOW_DEBUG
321       &       tempAdvScheme.NE.ENUM_CENTERED_2ND .AND.            IF ( debugLevel .GE. debLevB )
322       &       tempAdvScheme.NE.ENUM_UPWIND_3RD .AND.       &     CALL DEBUG_CALL('GAD_ADVECTION',myThid)
323       &       tempAdvScheme.NE.ENUM_CENTERED_4TH ) THEN  #endif
324            CALL GAD_ADVECTION(bi,bj,tempAdvScheme,GAD_TEMPERATURE,            CALL GAD_ADVECTION(
325       U                      theta,gT,       I             tempImplVertAdv, tempAdvScheme, tempVertAdvScheme,
326       I                      myTime,myIter,myThid)       I             GAD_TEMPERATURE,
327           ENDIF       I             uVel, vVel, wVel, theta,
328           IF (saltStepping .AND.       O             gT,
329       &       saltAdvScheme.NE.ENUM_CENTERED_2ND .AND.       I             bi,bj,myTime,myIter,myThid)
      &       saltAdvScheme.NE.ENUM_UPWIND_3RD .AND.  
      &       saltAdvScheme.NE.ENUM_CENTERED_4TH ) THEN  
           CALL GAD_ADVECTION(bi,bj,saltAdvScheme,GAD_SALINITY,  
      U                      salt,gS,  
      I                      myTime,myIter,myThid)  
          ENDIF  
330          ENDIF          ENDIF
331            IF (saltMultiDimAdvec) THEN
332    #ifdef ALLOW_DEBUG
333              IF ( debugLevel .GE. debLevB )
334         &     CALL DEBUG_CALL('GAD_ADVECTION',myThid)
335    #endif
336              CALL GAD_ADVECTION(
337         I             saltImplVertAdv, saltAdvScheme, saltVertAdvScheme,
338         I             GAD_SALINITY,
339         I             uVel, vVel, wVel, salt,
340         O             gS,
341         I             bi,bj,myTime,myIter,myThid)
342            ENDIF
343    C Since passive tracers are configurable separately from T,S we
344    C call the multi-dimensional method for PTRACERS regardless
345    C of whether multiDimAdvection is set or not.
346    #ifdef ALLOW_PTRACERS
347            IF ( usePTRACERS ) THEN
348    #ifdef ALLOW_DEBUG
349              IF ( debugLevel .GE. debLevB )
350         &     CALL DEBUG_CALL('PTRACERS_ADVECTION',myThid)
351    #endif
352             CALL PTRACERS_ADVECTION( bi,bj,myIter,myTime,myThid )
353            ENDIF
354    #endif /* ALLOW_PTRACERS */
355  #endif /* DISABLE_MULTIDIM_ADVECTION */  #endif /* DISABLE_MULTIDIM_ADVECTION */
356    
357    #ifdef ALLOW_DEBUG
358            IF ( debugLevel .GE. debLevB )
359         &    CALL DEBUG_MSG('ENTERING DOWNWARD K LOOP',myThid)
360    #endif
361    
362  C--     Start of thermodynamics loop  C--     Start of thermodynamics loop
363          DO k=Nr,1,-1          DO k=Nr,1,-1
364  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
365  C? Patrick Is this formula correct?  C? Patrick Is this formula correct?
366  cph Yes, but I rewrote it.  cph Yes, but I rewrote it.
367  cph Also, the KappaR? need the index and subscript k!  cph Also, the KappaR? need the index and subscript k!
368           kkey = (ikey-1)*Nr + k           kkey = (itdkey-1)*Nr + k
369  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
370    
371  C--       km1    Points to level above k (=k-1)  C--       km1    Points to level above k (=k-1)
# Line 500  C--       kDown  Cycles through 2,1 to p Line 381  C--       kDown  Cycles through 2,1 to p
381            jMin = 1-OLy            jMin = 1-OLy
382            jMax = sNy+OLy            jMax = sNy+OLy
383    
384              kp1Msk=1.
385              IF (k.EQ.Nr) kp1Msk=0.
386               DO j=1-Oly,sNy+Oly
387                DO i=1-Olx,sNx+Olx
388                 rTransKp1(i,j) = kp1Msk*rTrans(i,j)
389                ENDDO
390               ENDDO
391    #ifdef ALLOW_AUTODIFF_TAMC
392    CADJ STORE rTransKp1(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
393    #endif
394    
395  C--       Get temporary terms used by tendency routines  C--       Get temporary terms used by tendency routines
396            CALL CALC_COMMON_FACTORS (            CALL CALC_COMMON_FACTORS (
397       I         bi,bj,iMin,iMax,jMin,jMax,k,       I         bi,bj,iMin,iMax,jMin,jMax,k,
398       O         xA,yA,uTrans,vTrans,rTrans,maskUp,       O         xA,yA,uTrans,vTrans,rTrans,maskUp,
399       I         myThid)       I         myThid)
400    
401  #ifdef ALLOW_AUTODIFF_TAMC            IF (k.EQ.1) THEN
402  CADJ STORE KappaRT(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte  C- Surface interface :
403  CADJ STORE KappaRS(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte             DO j=1-Oly,sNy+Oly
404                DO i=1-Olx,sNx+Olx
405                 rTrans(i,j) = 0.
406                ENDDO
407               ENDDO
408              ELSE
409    C- Interior interface :
410               DO j=1-Oly,sNy+Oly
411                DO i=1-Olx,sNx+Olx
412                 rTrans(i,j) = rTrans(i,j)*maskC(i,j,k-1,bi,bj)
413                ENDDO
414               ENDDO
415              ENDIF
416    
417    #ifdef ALLOW_GMREDI
418    
419    C--   Residual transp = Bolus transp + Eulerian transp
420              IF (useGMRedi) THEN
421                CALL GMREDI_CALC_UVFLOW(
422         &            uTrans, vTrans, bi, bj, k, myThid)
423                IF (K.GE.2) CALL GMREDI_CALC_WFLOW(
424         &                    rTrans, bi, bj, k, myThid)
425              ENDIF
426    
427    #ifdef ALLOW_AUTODIFF_TAMC
428    CADJ STORE rTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
429    #ifdef GM_BOLUS_ADVEC
430    CADJ STORE uTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
431    CADJ STORE vTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
432    #endif
433  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
434    
435    #endif /* ALLOW_GMREDI */
436    
437  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL
438  C--      Calculate the total vertical diffusivity  C--      Calculate the total vertical diffusivity
439           CALL CALC_DIFFUSIVITY(           CALL CALC_DIFFUSIVITY(
# Line 518  C--      Calculate the total vertical di Line 441  C--      Calculate the total vertical di
441       I        maskUp,       I        maskUp,
442       O        KappaRT,KappaRS,       O        KappaRT,KappaRS,
443       I        myThid)       I        myThid)
444    # ifdef ALLOW_AUTODIFF_TAMC
445    CADJ STORE KappaRT(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte
446    CADJ STORE KappaRS(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte
447    # endif /* ALLOW_AUTODIFF_TAMC */
448  #endif  #endif
449    
450            iMin = 1-OLx+2            iMin = 1-OLx+2
# Line 530  C        and step forward storing result Line 457  C        and step forward storing result
457           IF ( tempStepping ) THEN           IF ( tempStepping ) THEN
458             CALL CALC_GT(             CALL CALC_GT(
459       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
460       I         xA,yA,uTrans,vTrans,rTrans,maskUp,       I         xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,
461       I         KappaRT,       I         KappaRT,
462       U         fVerT,       U         fVerT,
463       I         myTime,myIter,myThid)       I         myTime,myIter,myThid)
# Line 539  C        and step forward storing result Line 466  C        and step forward storing result
466       I         theta, gT,       I         theta, gT,
467       I         myIter, myThid)       I         myIter, myThid)
468           ENDIF           ENDIF
469    
470           IF ( saltStepping ) THEN           IF ( saltStepping ) THEN
471             CALL CALC_GS(             CALL CALC_GS(
472       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
473       I         xA,yA,uTrans,vTrans,rTrans,maskUp,       I         xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,
474       I         KappaRS,       I         KappaRS,
475       U         fVerS,       U         fVerS,
476       I         myTime,myIter,myThid)       I         myTime,myIter,myThid)
# Line 552  C        and step forward storing result Line 480  C        and step forward storing result
480       I         myIter, myThid)       I         myIter, myThid)
481           ENDIF           ENDIF
482  #ifdef ALLOW_PASSIVE_TRACER  #ifdef ALLOW_PASSIVE_TRACER
483    ceh3 needs an IF ( usePASSIVE_TRACER ) THEN
484           IF ( tr1Stepping ) THEN           IF ( tr1Stepping ) THEN
485             CALL CALC_GTR1(             CALL CALC_GTR1(
486       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
487       I         xA,yA,uTrans,vTrans,rTrans,maskUp,       I         xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,
488       I         KappaRT,       I         KappaRT,
489       U         fVerTr1,       U         fVerTr1,
490       I         myTime,myIter,myThid)       I         myTime,myIter,myThid)
# Line 565  C        and step forward storing result Line 494  C        and step forward storing result
494       I         myIter,myThid)       I         myIter,myThid)
495           ENDIF           ENDIF
496  #endif  #endif
497    #ifdef ALLOW_PTRACERS
498             IF ( usePTRACERS ) THEN
499               CALL PTRACERS_INTEGRATE(
500         I         bi,bj,k,
501         I         xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,
502         X         fVerP, KappaRS,
503         I         myIter,myTime,myThid)
504             ENDIF
505    #endif /* ALLOW_PTRACERS */
506    
507  #ifdef   ALLOW_OBCS  #ifdef   ALLOW_OBCS
508  C--      Apply open boundary conditions  C--      Apply open boundary conditions
# Line 574  C--      Apply open boundary conditions Line 512  C--      Apply open boundary conditions
512  #endif   /* ALLOW_OBCS */  #endif   /* ALLOW_OBCS */
513    
514  C--      Freeze water  C--      Freeze water
515           IF (allowFreezing) THEN  C  this bit of code is left here for backward compatibility.
516    C  freezing at surface level has been moved to FORWARD_STEP
517             IF ( useOldFreezing .AND. .NOT. useSEAICE
518         &       .AND. .NOT.(useThSIce.AND.k.EQ.1) ) THEN
519  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
520  CADJ STORE gT(:,:,k,bi,bj) = comlev1_bibj_k  CADJ STORE gT(:,:,k,bi,bj) = comlev1_bibj_k
521  CADJ &   , key = kkey, byte = isbyte  CADJ &   , key = kkey, byte = isbyte
522  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
523              CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )              CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )
524           END IF           ENDIF
525    
526  C--     end of thermodynamic k loop (Nr:1)  C--     end of thermodynamic k loop (Nr:1)
527          ENDDO          ENDDO
528    
529    
530    C--     Implicit vertical advection & diffusion
531    #ifdef INCLUDE_IMPLVERTADV_CODE
532            IF ( tempImplVertAdv ) THEN
533              CALL GAD_IMPLICIT_R(
534         I         tempImplVertAdv, tempAdvScheme, GAD_TEMPERATURE,
535         I         kappaRT, wVel, theta,
536         U         gT,
537         I         bi, bj, myTime, myIter, myThid )
538            ELSEIF ( tempStepping .AND. implicitDiffusion ) THEN
539    #else /* INCLUDE_IMPLVERTADV_CODE */
540            IF     ( tempStepping .AND. implicitDiffusion ) THEN
541    #endif /* INCLUDE_IMPLVERTADV_CODE */
542  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
543  C? Patrick? What about this one?  CADJ STORE KappaRT(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
544  cph Keys iikey and idkey dont seem to be needed  CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
 cph since storing occurs on different tape for each  
 cph impldiff call anyways.  
 cph Thus, common block comlev1_impl isnt needed either.  
 cph Storing below needed in the case useGMREDI.  
         iikey = (ikey-1)*maximpl  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 C--     Implicit diffusion  
         IF (implicitDiffusion) THEN  
   
          IF (tempStepping) THEN  
 #ifdef ALLOW_AUTODIFF_TAMC  
             idkey = iikey + 1  
 CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  
545  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
546              CALL IMPLDIFF(            CALL IMPLDIFF(
547       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
548       I         deltaTtracer, KappaRT, recip_HFacC,       I         deltaTtracer, KappaRT, recip_HFacC,
549       U         gT,       U         gT,
550       I         myThid )       I         myThid )
551           ENDIF          ENDIF
552    
553           IF (saltStepping) THEN  #ifdef INCLUDE_IMPLVERTADV_CODE
554            IF ( saltImplVertAdv ) THEN
555              CALL GAD_IMPLICIT_R(
556         I         saltImplVertAdv, saltAdvScheme, GAD_SALINITY,
557         I         kappaRS, wVel, salt,
558         U         gS,
559         I         bi, bj, myTime, myIter, myThid )
560            ELSEIF ( saltStepping .AND. implicitDiffusion ) THEN
561    #else /* INCLUDE_IMPLVERTADV_CODE */
562            IF     ( saltStepping .AND. implicitDiffusion ) THEN
563    #endif /* INCLUDE_IMPLVERTADV_CODE */
564  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
565           idkey = iikey + 2  CADJ STORE KappaRS(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
566  CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
567  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
568              CALL IMPLDIFF(            CALL IMPLDIFF(
569       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
570       I         deltaTtracer, KappaRS, recip_HFacC,       I         deltaTtracer, KappaRS, recip_HFacC,
571       U         gS,       U         gS,
572       I         myThid )       I         myThid )
573           ENDIF          ENDIF
574    
575  #ifdef ALLOW_PASSIVE_TRACER  #ifdef ALLOW_PASSIVE_TRACER
576           IF (tr1Stepping) THEN          IF ( tr1Stepping .AND. implicitDiffusion ) THEN
577  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
578  CADJ STORE gTr1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE gTr1(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
579  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
580            CALL IMPLDIFF(            CALL IMPLDIFF(
581       I      bi, bj, iMin, iMax, jMin, jMax,       I      bi, bj, iMin, iMax, jMin, jMax,
582       I      deltaTtracer, KappaRT, recip_HFacC,       I      deltaTtracer, KappaRT, recip_HFacC,
583       U      gTr1,       U      gTr1,
584       I      myThid )       I      myThid )
585           ENDIF          ENDIF
586  #endif  #endif
587    
588    #ifdef ALLOW_PTRACERS
589    c #ifdef INCLUDE_IMPLVERTADV_CODE
590    c       IF ( usePTRACERS .AND. ptracerImplVertAdv ) THEN
591    c       ELSEIF ( usePTRACERS .AND. implicitDiffusion ) THEN
592    c #else
593            IF     ( usePTRACERS .AND. implicitDiffusion ) THEN
594    C--     Vertical diffusion (implicit) for passive tracers
595               CALL PTRACERS_IMPLDIFF( bi,bj,KappaRS,myThid )
596            ENDIF
597    #endif /* ALLOW_PTRACERS */
598    
599  #ifdef   ALLOW_OBCS  #ifdef   ALLOW_OBCS
600  C--      Apply open boundary conditions  C--      Apply open boundary conditions
601           IF (useOBCS) THEN          IF ( ( implicitDiffusion
602         &    .OR. tempImplVertAdv
603         &    .OR. saltImplVertAdv
604         &       ) .AND. useOBCS     ) THEN
605             DO K=1,Nr             DO K=1,Nr
606               CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )               CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )
607             ENDDO             ENDDO
608           END IF          ENDIF
609  #endif   /* ALLOW_OBCS */  #endif   /* ALLOW_OBCS */
610    
611  C--     End If implicitDiffusion  #ifdef ALLOW_TIMEAVE
612    #ifndef HRCUBE
613            IF (taveFreq.GT.0. .AND. ivdc_kappa.NE.0.) THEN
614              CALL TIMEAVE_CUMULATE(ConvectCountTave, IVDConvCount,
615         I                           Nr, deltaTclock, bi, bj, myThid)
616          ENDIF          ENDIF
617            useVariableK = useKPP .OR. useGMredi .OR. ivdc_kappa.NE.0.
618            IF (taveFreq.GT.0. .AND. useVariableK ) THEN
619             IF (implicitDiffusion) THEN
620              CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, gT, kappaRT,
621         I                        Nr, 3, deltaTclock, bi, bj, myThid)
622             ELSE
623              CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, theta, kappaRT,
624         I                        Nr, 3, deltaTclock, bi, bj, myThid)
625             ENDIF
626            ENDIF
627    #endif /* ndef HRCUBE */
628    #endif /* ALLOW_TIMEAVE */
629    
630    #endif /* SINGLE_LAYER_MODE */
631    
632  Ccs-  C--   end bi,bj loops.
633         ENDDO         ENDDO
634        ENDDO        ENDDO
635    
636  #ifdef ALLOW_AIM  #ifdef ALLOW_DEBUG
637        IF ( useAIM ) THEN        If (debugMode) THEN
638         CALL AIM_AIM2DYN_EXCHANGES( myTime, myThid )         CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (THERMODYNAMICS)',myThid)
639           CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (THERMODYNAMICS)',myThid)
640           CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (THERMODYNAMICS)',myThid)
641           CALL DEBUG_STATS_RL(Nr,theta,'Theta (THERMODYNAMICS)',myThid)
642           CALL DEBUG_STATS_RL(Nr,salt,'Salt (THERMODYNAMICS)',myThid)
643           CALL DEBUG_STATS_RL(Nr,Gt,'Gt (THERMODYNAMICS)',myThid)
644           CALL DEBUG_STATS_RL(Nr,Gs,'Gs (THERMODYNAMICS)',myThid)
645           CALL DEBUG_STATS_RL(Nr,GtNm1,'GtNm1 (THERMODYNAMICS)',myThid)
646           CALL DEBUG_STATS_RL(Nr,GsNm1,'GsNm1 (THERMODYNAMICS)',myThid)
647    #ifdef ALLOW_PTRACERS
648           IF ( usePTRACERS ) THEN
649             CALL PTRACERS_DEBUG(myThid)
650           ENDIF
651    #endif /* ALLOW_PTRACERS */
652        ENDIF        ENDIF
653         _EXCH_XYZ_R8(gT,myThid)  #endif
654         _EXCH_XYZ_R8(gS,myThid)  
655  #else  #ifdef ALLOW_DEBUG
656        IF (staggerTimeStep.AND.useCubedSphereExchange) THEN           IF ( debugLevel .GE. debLevB )
657         _EXCH_XYZ_R8(gT,myThid)       &    CALL DEBUG_LEAVE('THERMODYNAMICS',myThid)
658         _EXCH_XYZ_R8(gS,myThid)  #endif
       ENDIF  
 #endif /* ALLOW_AIM */  
659    
660        RETURN        RETURN
661        END        END
Legend:
Removed from v.1.14  
changed lines
  Added in v.1.72
  ViewVC Help
Powered by ViewVC 1.1.22