/[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.36 by heimbach, Tue Jan 21 19:19:45 2003 UTC revision 1.73 by jmc, Wed Jul 7 01:32:11 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  #ifdef ALLOW_AUTODIFF_TAMC
11  # ifdef ALLOW_GMREDI  # ifdef ALLOW_GMREDI
12  #  include "GMREDI_OPTIONS.h"  #  include "GMREDI_OPTIONS.h"
# Line 9  C $Name$ Line 14  C $Name$
14  # ifdef ALLOW_KPP  # ifdef ALLOW_KPP
15  #  include "KPP_OPTIONS.h"  #  include "KPP_OPTIONS.h"
16  # endif  # endif
 cswdice --- add ----  
 #ifdef ALLOW_THERM_SEAICE  
 #include "ICE.h"  
 #endif  
 cswdice ------  
17  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
18    
19  CBOP  CBOP
# Line 84  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"
# Line 95  C     == Global variables === Line 102  C     == Global variables ===
102  # ifdef ALLOW_GMREDI  # ifdef ALLOW_GMREDI
103  #  include "GMREDI.h"  #  include "GMREDI.h"
104  # endif  # endif
105    # ifdef ALLOW_EBM
106    #  include "EBM.h"
107    # endif
108  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
 #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 117  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 143  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
# Line 173  C--   dummy statement to end declaration Line 181  C--   dummy statement to end declaration
181        itdkey = 1        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  
         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 201  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 221  CHPF$&                  ) Line 211  CHPF$&                  )
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            rhoKM1 (i,j)   = 0. _d 0  #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 239  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
            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  
            ConvectCount(i,j,k) = 0.  
249             KappaRT(i,j,k)    = 0. _d 0             KappaRT(i,j,k)    = 0. _d 0
250             KappaRS(i,j,k)    = 0. _d 0             KappaRS(i,j,k)    = 0. _d 0
251  #ifdef ALLOW_AUTODIFF_TAMC  C- tracer tendency needs to be set to zero (moved here from gad_calc_rhs):
 cph all the following init. are necessary for TAF  
 cph although some of these are re-initialised later.  
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
254  # ifdef ALLOW_PASSIVE_TRACER  # ifdef ALLOW_PASSIVE_TRACER
255    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  # ifdef ALLOW_GMREDI  # ifdef ALLOW_PTRACERS
259             Kwx(i,j,k,bi,bj)  = 0. _d 0  ceh3 this should have an   IF ( usePTRACERS ) THEN
260             Kwy(i,j,k,bi,bj)  = 0. _d 0             DO iTracer=1,PTRACERS_numInUse
261             Kwz(i,j,k,bi,bj)  = 0. _d 0              gPTr(i,j,k,bi,bj,itracer) = 0. _d 0
262  #  ifdef GM_NON_UNITY_DIAGONAL             ENDDO
263             Kux(i,j,k,bi,bj)  = 0. _d 0  # endif
            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 */  
264            ENDDO            ENDDO
265           ENDDO           ENDDO
266          ENDDO          ENDDO
267    
268          iMin = 1-OLx  c       iMin = 1-OLx
269          iMax = sNx+OLx  c       iMax = sNx+OLx
270          jMin = 1-OLy  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=itdkey, byte=isbyte  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
275  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, 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=itdkey, byte=isbyte  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
280  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, 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 = (itdkey-1)*Nr + k  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 C--       Integrate continuity vertically for vertical velocity  
 c         CALL INTEGRATE_FOR_W(  
 c    I                         bi, bj, k, uVel, vVel,  
 c    O                         wVel,  
 c    I                         myThid )  
   
 #ifdef    ALLOW_OBCS  
 #ifdef    ALLOW_NONHYDROSTATIC  
 C--       Apply OBC to W if in N-H mode  
 c         IF (useOBCS.AND.nonHydrostatic) THEN  
 c           CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )  
 c         ENDIF  
 #endif    /* ALLOW_NONHYDROSTATIC */  
 #endif    /* ALLOW_OBCS */  
   
 C--       Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h  
 C--       MOST of THERMODYNAMICS will be disabled  
 #ifndef SINGLE_LAYER_MODE  
   
 C--       Calculate gradients of potential density for isoneutral  
 C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)  
 c         IF ( k.GT.1 .AND. (useGMRedi.OR.ivdc_kappa.NE.0.) ) THEN  
           IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.) ) THEN  
 #ifdef ALLOW_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  
 CADJ STORE pressure(:,:,k,bi,bj) =  
 CADJ &     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  
 CADJ STORE pressure(:,:,k-1,bi,bj) =  
 CADJ &     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  
             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  
             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  
   
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=itdkey, byte=isbyte  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
 CADJ STORE pressure (:,:,:,bi,bj) =  
 CADJ &     comlev1_bibj, key=itdkey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 #ifdef  ALLOW_OBCS  
 C--     Calculate future values on open boundaries  
         IF (useOBCS) THEN  
           CALL OBCS_CALC( bi, bj, myTime+deltaT, myIter+1,  
      I            uVel, vVel, wVel, theta, salt,  
      I            myThid )  
         ENDIF  
 #endif  /* ALLOW_OBCS */  
   
   
 c********************************************  
 cswdice --- add ---  
 #ifdef ALLOW_THERM_SEAICE  
 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  
   
 cswdice --- end add ---  
   
 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 )  
 cswdice --- add ----  
 #endif  
 cswdice --- end add ---  
 c******************************************  
   
   
   
   
 #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  
287  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
288    
289  C--     Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h  C--     Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h
290  C--     MOST of THERMODYNAMICS will be disabled  C--     MOST of THERMODYNAMICS will be disabled
291  #ifndef SINGLE_LAYER_MODE  #ifndef SINGLE_LAYER_MODE
292    
 #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  
           CALL GMREDI_CALC_TENSOR(  
      I             bi, bj, iMin, iMax, jMin, jMax,  
      I             sigmaX, sigmaY, sigmaR,  
      I             myThid )  
293  #ifdef ALLOW_AUTODIFF_TAMC  #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  
           CALL KPP_CALC(  
      I                  bi, bj, myTime, myThid )  
 #ifdef ALLOW_AUTODIFF_TAMC  
         ELSE  
           CALL KPP_CALC_DUMMY(  
      I                  bi, bj, myTime, myThid )  
 #endif /* ALLOW_AUTODIFF_TAMC */  
         ENDIF  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE KPPghat   (:,:,:,bi,bj)  
 CADJ &   , KPPdiffKzT(:,:,:,bi,bj)  
 CADJ &   , KPPdiffKzS(:,:,:,bi,bj)  
 CADJ &   , KPPfrac   (:,:  ,bi,bj)  
 CADJ &                 = comlev1_bibj, key=itdkey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 #endif  /* ALLOW_KPP */  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE KappaRT(:,:,:)     = comlev1_bibj, key=itdkey, byte=isbyte  
 CADJ STORE KappaRS(:,:,:)     = comlev1_bibj, key=itdkey, byte=isbyte  
294  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
295  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
296  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
# Line 503  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_ Line 298  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_
298  #ifdef ALLOW_PASSIVE_TRACER  #ifdef ALLOW_PASSIVE_TRACER
299  CADJ STORE tr1  (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE tr1  (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
300  #endif  #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
306  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
307    
 #ifdef ALLOW_AIM  
 C       AIM - atmospheric intermediate model, physics package code.  
         IF ( useAIM ) THEN  
          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 */  
   
 #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 533  C recomputation. It *is* differentiable, Line 317  C recomputation. It *is* differentiable,
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 (tempMultiDimAdvec) THEN          IF (tempMultiDimAdvec) THEN
320            CALL GAD_ADVECTION(bi,bj,tempAdvScheme,GAD_TEMPERATURE,  #ifdef ALLOW_DEBUG
321       U                      theta,gT,            IF ( debugLevel .GE. debLevB )
322       I                      myTime,myIter,myThid)       &     CALL DEBUG_CALL('GAD_ADVECTION',myThid)
323    #endif
324              CALL GAD_ADVECTION(
325         I             tempImplVertAdv, tempAdvScheme, tempVertAdvScheme,
326         I             GAD_TEMPERATURE,
327         I             uVel, vVel, wVel, theta,
328         O             gT,
329         I             bi,bj,myTime,myIter,myThid)
330          ENDIF          ENDIF
331          IF (saltMultiDimAdvec) THEN          IF (saltMultiDimAdvec) THEN
332            CALL GAD_ADVECTION(bi,bj,saltAdvScheme,GAD_SALINITY,  #ifdef ALLOW_DEBUG
333       U                      salt,gS,            IF ( debugLevel .GE. debLevB )
334       I                      myTime,myIter,myThid)       &     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          ENDIF
343  C Since passive tracers are configurable separately from T,S we  C Since passive tracers are configurable separately from T,S we
344  C call the multi-dimensional method for PTRACERS regardless  C call the multi-dimensional method for PTRACERS regardless
345  C of whether multiDimAdvection is set or not.  C of whether multiDimAdvection is set or not.
346  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
347          IF ( usePTRACERS ) THEN          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 )           CALL PTRACERS_ADVECTION( bi,bj,myIter,myTime,myThid )
353          ENDIF          ENDIF
354  #endif /* ALLOW_PTRACERS */  #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
# Line 574  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              IF (k.EQ.1) THEN
402    C- Surface interface :
403               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  #ifdef ALLOW_GMREDI
418    
419  C--   Residual transp = Bolus transp + Eulerian transp  C--   Residual transp = Bolus transp + Eulerian transp
# Line 590  C--   Residual transp = Bolus transp + E Line 424  C--   Residual transp = Bolus transp + E
424       &                    rTrans, bi, bj, k, myThid)       &                    rTrans, bi, bj, k, myThid)
425            ENDIF            ENDIF
426    
427  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
428    CADJ STORE rTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
429  #ifdef GM_BOLUS_ADVEC  #ifdef GM_BOLUS_ADVEC
430  CADJ STORE uTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE uTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
431  CADJ STORE vTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE vTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
 CADJ STORE rTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte  
432  #endif  #endif
433  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
434    
435  #endif /* ALLOW_GMREDI */  #endif /* ALLOW_GMREDI */
436    
 #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 */  
   
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 612  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 624  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 633  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  cswdice ---- add ---  
 #ifdef ALLOW_THERM_SEAICE  
        if (k.eq.1) then  
         call ICE_FREEZE(bi, bj, iMin, iMax, jMin, jMax, myThid )  
        endif  
 #endif  
 cswdice -- end add ---  
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 653  cswdice -- end add --- Line 480  cswdice -- end add ---
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 668  cswdice -- end add --- Line 496  cswdice -- end add ---
496  #endif  #endif
497  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
498           IF ( usePTRACERS ) THEN           IF ( usePTRACERS ) THEN
499             CALL PTRACERS_INTEGERATE(             CALL PTRACERS_INTEGRATE(
500       I         bi,bj,k,       I         bi,bj,k,
501       I         xA,yA,uTrans,vTrans,rTrans,maskUp,       I         xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,
502       X         KappaRS,       X         fVerP, KappaRS,
503       I         myIter,myTime,myThid)       I         myIter,myTime,myThid)
504           ENDIF           ENDIF
505  #endif /* ALLOW_PTRACERS */  #endif /* ALLOW_PTRACERS */
# Line 684  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 .AND. .NOT. useSEAICE ) 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    
 cswdice -- add ---  
 #ifdef ALLOW_THERM_SEAICE  
 c timeaveraging for ice model values  
            CALL ICE_AVE(bi,bj,iMin,iMax,jMin,jMax,myThid )  
 #endif  
 cswdice --- end add ---  
   
   
   
   
 C--     Implicit diffusion  
         IF (implicitDiffusion) THEN  
529    
530           IF (tempStepping) THEN  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    CADJ STORE KappaRT(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
544  CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte  CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, 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    CADJ STORE KappaRS(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
566  CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, 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=itdkey, byte=isbyte  CADJ STORE gTr1(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
579  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
# Line 740  CADJ STORE gTr1(:,:,:,bi,bj) = comlev1_b Line 582  CADJ STORE gTr1(:,:,:,bi,bj) = comlev1_b
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  #ifdef ALLOW_PTRACERS
589  C Vertical diffusion (implicit) for passive tracers  c #ifdef INCLUDE_IMPLVERTADV_CODE
590           IF ( usePTRACERS ) THEN  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 )             CALL PTRACERS_IMPLDIFF( bi,bj,KappaRS,myThid )
596           ENDIF          ENDIF
597  #endif /* ALLOW_PTRACERS */  #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            IF ( taveFreq.GT. 0. _d 0 .AND.
613         &       buoyancyRelation(1:7) .EQ. 'OCEANIC' ) THEN
614              CALL TIMEAVE_SURF_FLUX( bi, bj, myTime, myIter, myThid)
615            ENDIF
616    #ifndef HRCUBE
617            IF (taveFreq.GT.0. .AND. ivdc_kappa.NE.0.) THEN
618              CALL TIMEAVE_CUMULATE(ConvectCountTave, IVDConvCount,
619         I                           Nr, deltaTclock, bi, bj, myThid)
620            ENDIF
621            useVariableK = useKPP .OR. useGMredi .OR. ivdc_kappa.NE.0.
622            IF (taveFreq.GT.0. .AND. useVariableK ) THEN
623             IF (implicitDiffusion) THEN
624              CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, gT, kappaRT,
625         I                        Nr, 3, deltaTclock, bi, bj, myThid)
626             ELSE
627              CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, theta, kappaRT,
628         I                        Nr, 3, deltaTclock, bi, bj, myThid)
629             ENDIF
630          ENDIF          ENDIF
631    #endif /* ndef HRCUBE */
632    #endif /* ALLOW_TIMEAVE */
633    
634  #endif /* SINGLE_LAYER_MODE */  #endif /* SINGLE_LAYER_MODE */
635    
636  Ccs-  C--   end bi,bj loops.
637         ENDDO         ENDDO
638        ENDDO        ENDDO
639    
640  #ifdef ALLOW_AIM  #ifdef ALLOW_DEBUG
 c     IF ( useAIM ) THEN  
 c      CALL AIM_AIM2DYN_EXCHANGES( myTime, myThid )  
 c     ENDIF  
 #endif /* ALLOW_AIM */  
 c     IF ( staggerTimeStep ) THEN  
 c      IF ( useAIM .OR. useCubedSphereExchange ) THEN  
 c        IF (tempStepping) _EXCH_XYZ_R8(gT,myThid)  
 c        IF (saltStepping) _EXCH_XYZ_R8(gS,myThid)  
 c      ELSEIF ( useGMRedi .AND. Oly.LT.4 ) THEN  
 cc       .AND. GM_AdvForm .AND. .NOT.GM_AdvSeparate ) THEN  
 c        IF (tempMultiDimAdvec) _EXCH_XYZ_R8(gT,myThid)  
 c        IF (saltMultiDimAdvec) _EXCH_XYZ_R8(gS,myThid)  
 c      ENDIF    
 c     ENDIF    
   
 #ifndef DISABLE_DEBUGMODE  
641        If (debugMode) THEN        If (debugMode) THEN
642         CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (THERMODYNAMICS)',myThid)
643         CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (THERMODYNAMICS)',myThid)
# Line 803  c     ENDIF Line 656  c     ENDIF
656        ENDIF        ENDIF
657  #endif  #endif
658    
659    #ifdef ALLOW_DEBUG
660             IF ( debugLevel .GE. debLevB )
661         &    CALL DEBUG_LEAVE('THERMODYNAMICS',myThid)
662    #endif
663    
664        RETURN        RETURN
665        END        END
Legend:
Removed from v.1.36  
changed lines
  Added in v.1.73
  ViewVC Help
Powered by ViewVC 1.1.22