/[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.24 by jmc, Wed Jul 3 20:22:39 2002 UTC revision 1.91 by spk, Wed Apr 20 15:52:22 2005 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 75  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_OFFLINE
84    #include "OFFLINE.h"
85    #endif
86    #ifdef ALLOW_PTRACERS
87    #include "PTRACERS_SIZE.h"
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 */
109  #ifdef ALLOW_TIMEAVE  
110  #include "TIMEAVE_STATV.h"  #ifdef ALLOW_MATRIX
111  #endif  # include "MATRIX.h"
112    #endif /* ALLOW_MATRIX */
113    
114  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
115  C     == Routine arguments ==  C     == Routine arguments ==
# Line 103  C     myThid - Thread number for this in Line 120  C     myThid - Thread number for this in
120        INTEGER myIter        INTEGER myIter
121        INTEGER myThid        INTEGER myThid
122    
123    #ifdef ALLOW_GENERIC_ADVDIFF
124  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
125  C     == Local variables  C     == Local variables
126  C     xA, yA                 - Per block temporaries holding face areas  C     xA, yA                 - Per block temporaries holding face areas
# Line 111  C                              transport Line 129  C                              transport
129  C                              o uTrans: Zonal transport  C                              o uTrans: Zonal transport
130  C                              o vTrans: Meridional transport  C                              o vTrans: Meridional transport
131  C                              o rTrans: Vertical transport  C                              o rTrans: Vertical transport
132    C     rTransKp1                o vertical volume transp. at interface k+1
133  C     maskUp                   o maskUp: land/water mask for W points  C     maskUp                   o maskUp: land/water mask for W points
134  C     fVer[STUV]               o fVer: Vertical flux term - note fVer  C     fVer[STUV]               o fVer: Vertical flux term - note fVer
135  C                                      is "pipelined" in the vertical  C                                      is "pipelined" in the vertical
136  C                                      so we need an fVer for each  C                                      so we need an fVer for each
137  C                                      variable.  C                                      variable.
138  C     rhoK, rhoKM1   - Density at current level, and level above  C     kappaRT,       - Total diffusion in vertical at level k, for T and S
139  C     phiHyd         - Hydrostatic part of the potential phiHydi.  C     kappaRS          (background + spatially varying, isopycnal term).
140  C                      In z coords phiHydiHyd is the hydrostatic  C     kappaRTr       - Total diffusion in vertical at level k,
141  C                      Potential (=pressure/rho0) anomaly  C                      for each passive Tracer
142  C                      In p coords phiHydiHyd is the geopotential  C     kappaRk        - Total diffusion in vertical, all levels, 1 tracer  
143  C                      surface height anomaly.  C     useVariableK   = T when vertical diffusion is not constant
 C     phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean)  
 C     phiSurfY             or geopotentiel (atmos) in X and Y direction  
 C     KappaRT,       - Total diffusion in vertical for T and S.  
 C     KappaRS          (background + spatially varying, isopycnal term).  
144  C     iMin, iMax     - Ranges and sub-block indices on which calculations  C     iMin, iMax     - Ranges and sub-block indices on which calculations
145  C     jMin, jMax       are applied.  C     jMin, jMax       are applied.
146  C     bi, bj  C     bi, bj
# Line 137  C                      index into fVerTe Line 152  C                      index into fVerTe
152        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
153        _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
154        _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
155          _RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
156        _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
157        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
158        _RL fVerS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
159        _RL fVerTr1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL kappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
160        _RL phiHyd  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL kappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
161        _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  #ifdef ALLOW_PTRACERS
162        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fVerP   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2,PTRACERS_num)
163        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL kappaRTr(1-Olx:sNx+Olx,1-Oly:sNy+Oly,PTRACERS_num)
164        _RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  #endif
165        _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL kappaRk (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
166        _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL kp1Msk
167        _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        LOGICAL useVariableK
       _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)  
168        INTEGER iMin, iMax        INTEGER iMin, iMax
169        INTEGER jMin, jMax        INTEGER jMin, jMax
170        INTEGER bi, bj        INTEGER bi, bj
171        INTEGER i, j        INTEGER i, j
172        INTEGER k, km1, kup, kDown        INTEGER k, km1, kup, kDown
173          INTEGER iTracer, ip
174    
175  CEOP  CEOP
176    
177    #ifdef ALLOW_DEBUG
178             IF ( debugLevel .GE. debLevB )
179         &    CALL DEBUG_ENTER('THERMODYNAMICS',myThid)
180    #endif
181    
182  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
183  C--   dummy statement to end declaration part  C--   dummy statement to end declaration part
184        ikey = 1        ikey = 1
185          itdkey = 1
186  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
187    
 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  
   
   
188  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
189  C--   HPF directive to help TAMC  C--   HPF directive to help TAMC
190  CHPF$ INDEPENDENT  CHPF$ INDEPENDENT
# Line 194  CHPF$ INDEPENDENT Line 195  CHPF$ INDEPENDENT
195  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
196  C--    HPF directive to help TAMC  C--    HPF directive to help TAMC
197  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS
198  CHPF$&                  ,phiHyd,utrans,vtrans,xA,yA  CHPF$&                  ,utrans,vtrans,xA,yA
199  CHPF$&                  ,KappaRT,KappaRS  CHPF$&                  ,kappaRT,kappaRS
200  CHPF$&                  )  CHPF$&                  )
201  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
202    
# Line 209  CHPF$&                  ) Line 210  CHPF$&                  )
210            act3 = myThid - 1            act3 = myThid - 1
211            max3 = nTx*nTy            max3 = nTx*nTy
212            act4 = ikey_dynamics - 1            act4 = ikey_dynamics - 1
213            ikey = (act1 + 1) + act2*max1            itdkey = (act1 + 1) + act2*max1
214       &                      + act3*max1*max2       &                      + act3*max1*max2
215       &                      + act4*max1*max2*max3       &                      + act4*max1*max2*max3
216  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
217    
218  C--     Set up work arrays that need valid initial values  C--   Set up work arrays with valid (i.e. not NaN) values
219    C     These inital values do not alter the numerical results. They
220    C     just ensure that all memory references are to valid floating
221    C     point numbers. This prevents spurious hardware signals due to
222    C     uninitialised but inert locations.
223    
224          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
225           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
226              xA(i,j)        = 0. _d 0
227              yA(i,j)        = 0. _d 0
228              uTrans(i,j)    = 0. _d 0
229              vTrans(i,j)    = 0. _d 0
230            rTrans (i,j)   = 0. _d 0            rTrans (i,j)   = 0. _d 0
231              rTransKp1(i,j) = 0. _d 0
232            fVerT  (i,j,1) = 0. _d 0            fVerT  (i,j,1) = 0. _d 0
233            fVerT  (i,j,2) = 0. _d 0            fVerT  (i,j,2) = 0. _d 0
234            fVerS  (i,j,1) = 0. _d 0            fVerS  (i,j,1) = 0. _d 0
235            fVerS  (i,j,2) = 0. _d 0            fVerS  (i,j,2) = 0. _d 0
236            fVerTr1(i,j,1) = 0. _d 0            kappaRT(i,j)   = 0. _d 0
237            fVerTr1(i,j,2) = 0. _d 0            kappaRS(i,j)   = 0. _d 0
238            rhoKM1 (i,j)   = 0. _d 0  #ifdef ALLOW_PTRACERS
239              DO ip=1,PTRACERS_num
240               fVerP  (i,j,1,ip) = 0. _d 0
241               fVerP  (i,j,2,ip) = 0. _d 0
242               kappaRTr(i,j,ip)  = 0. _d 0
243              ENDDO
244    #endif
245           ENDDO           ENDDO
246          ENDDO          ENDDO
247    
# Line 232  C--     Set up work arrays that need val Line 249  C--     Set up work arrays that need val
249           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
250            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
251  C This is currently also used by IVDC and Diagnostics  C This is currently also used by IVDC and Diagnostics
252             phiHyd(i,j,k)  = 0. _d 0             kappaRk(i,j,k)    = 0. _d 0
253             sigmaX(i,j,k) = 0. _d 0  C- tracer tendency needs to be set to zero (moved here from gad_calc_rhs):
254             sigmaY(i,j,k) = 0. _d 0             gT(i,j,k,bi,bj)   = 0. _d 0
255             sigmaR(i,j,k) = 0. _d 0             gS(i,j,k,bi,bj)   = 0. _d 0
256             ConvectCount(i,j,k) = 0.  # ifdef ALLOW_PTRACERS
257             KappaRT(i,j,k) = 0. _d 0             DO iTracer=1,PTRACERS_numInUse
258             KappaRS(i,j,k) = 0. _d 0              gPTr(i,j,k,bi,bj,itracer) = 0. _d 0
259  #ifdef ALLOW_AUTODIFF_TAMC             ENDDO
260             gT(i,j,k,bi,bj) = 0. _d 0  # endif
            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_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  
 #endif /* ALLOW_GMREDI */  
 #endif  
261            ENDDO            ENDDO
262           ENDDO           ENDDO
263          ENDDO          ENDDO
264    
265          iMin = 1-OLx  c       iMin = 1-OLx
266          iMax = sNx+OLx  c       iMax = sNx+OLx
267          jMin = 1-OLy  c       jMin = 1-OLy
268          jMax = sNy+OLy  c       jMax = sNy+OLy
   
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  
 #ifdef ALLOW_KPP  
 CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  
 #endif  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 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--       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  
 #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  
   
 #endif /* SINGLE_LAYER_MODE */  
   
 C--     end of diagnostic k loop (Nr:1)  
         ENDDO  
269    
270  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
271  cph avoids recomputation of integrate_for_w  cph avoids recomputation of integrate_for_w
272  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE wvel (:,:,:,bi,bj) = 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--     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  
273  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
274    
275  C--     Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h  C--     Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h
276  C--     MOST of THERMODYNAMICS will be disabled  C--     MOST of THERMODYNAMICS will be disabled
277  #ifndef SINGLE_LAYER_MODE  #ifndef SINGLE_LAYER_MODE
278    
 #ifdef  ALLOW_GMREDI  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE sigmaX(:,:,k)        = comlev1_bibj_k, key=kkey, byte=isbyte  
 CADJ STORE sigmaY(:,:,k)        = comlev1_bibj_k, key=kkey, byte=isbyte  
 CADJ STORE sigmaR(:,:,k)        = comlev1_bibj_k, key=kkey, 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 )  
 #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=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  
   
 #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=ikey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 #endif  /* ALLOW_KPP */  
   
279  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
280  CADJ STORE KappaRT(:,:,:)     = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
281  CADJ STORE KappaRS(:,:,:)     = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
282  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
283  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
284  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  #ifdef ALLOW_PTRACERS
285  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  cph-- moved to forward_step to avoid key computation
286  #ifdef ALLOW_PASSIVE_TRACER  cphCADJ STORE ptracer(:,:,:,bi,bj,itracer) = comlev1_bibj,
287  CADJ STORE tr1  (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  cphCADJ &                              key=itdkey, byte=isbyte
288  #endif  #endif
289  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
290    
 #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 */  
   
291  #ifndef DISABLE_MULTIDIM_ADVECTION  #ifndef DISABLE_MULTIDIM_ADVECTION
292  C--     Some advection schemes are better calculated using a multi-dimensional  C--     Some advection schemes are better calculated using a multi-dimensional
293  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 479  C to be able to exclude this scheme to a Line 299  C to be able to exclude this scheme to a
299  C recomputation. It *is* differentiable, if you need it.  C recomputation. It *is* differentiable, if you need it.
300  C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to  C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to
301  C disable this section of code.  C disable this section of code.
302    #ifndef ALLOW_OFFLINE
303          IF (tempMultiDimAdvec) THEN          IF (tempMultiDimAdvec) THEN
304            CALL GAD_ADVECTION(bi,bj,tempAdvScheme,GAD_TEMPERATURE,  #ifdef ALLOW_DEBUG
305       U                      theta,gT,            IF ( debugLevel .GE. debLevB )
306       I                      myTime,myIter,myThid)       &     CALL DEBUG_CALL('GAD_ADVECTION',myThid)
307    #endif
308              CALL GAD_ADVECTION(
309         I             tempImplVertAdv, tempAdvScheme, tempVertAdvScheme,
310         I             GAD_TEMPERATURE,
311         I             uVel, vVel, wVel, theta,
312         O             gT,
313         I             bi,bj,myTime,myIter,myThid)
314          ENDIF          ENDIF
315    #endif
316    #ifndef ALLOW_OFFLINE
317          IF (saltMultiDimAdvec) THEN          IF (saltMultiDimAdvec) THEN
318            CALL GAD_ADVECTION(bi,bj,saltAdvScheme,GAD_SALINITY,  #ifdef ALLOW_DEBUG
319       U                      salt,gS,            IF ( debugLevel .GE. debLevB )
320       I                      myTime,myIter,myThid)       &     CALL DEBUG_CALL('GAD_ADVECTION',myThid)
321    #endif
322              CALL GAD_ADVECTION(
323         I             saltImplVertAdv, saltAdvScheme, saltVertAdvScheme,
324         I             GAD_SALINITY,
325         I             uVel, vVel, wVel, salt,
326         O             gS,
327         I             bi,bj,myTime,myIter,myThid)
328          ENDIF          ENDIF
329    #endif
330  C Since passive tracers are configurable separately from T,S we  C Since passive tracers are configurable separately from T,S we
331  C call the multi-dimensional method for PTRACERS regardless  C call the multi-dimensional method for PTRACERS regardless
332  C of whether multiDimAdvection is set or not.  C of whether multiDimAdvection is set or not.
333  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
334          IF ( usePTRACERS ) THEN          IF ( usePTRACERS ) THEN
335    #ifdef ALLOW_DEBUG
336              IF ( debugLevel .GE. debLevB )
337         &     CALL DEBUG_CALL('PTRACERS_ADVECTION',myThid)
338    #endif
339           CALL PTRACERS_ADVECTION( bi,bj,myIter,myTime,myThid )           CALL PTRACERS_ADVECTION( bi,bj,myIter,myTime,myThid )
340          ENDIF          ENDIF
341  #endif /* ALLOW_PTRACERS */  #endif /* ALLOW_PTRACERS */
342  #endif /* DISABLE_MULTIDIM_ADVECTION */  #endif /* DISABLE_MULTIDIM_ADVECTION */
343    
344    #ifdef ALLOW_DEBUG
345            IF ( debugLevel .GE. debLevB )
346         &    CALL DEBUG_MSG('ENTERING DOWNWARD K LOOP',myThid)
347    #endif
348    
349  C--     Start of thermodynamics loop  C--     Start of thermodynamics loop
350          DO k=Nr,1,-1          DO k=Nr,1,-1
351  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
352  C? Patrick Is this formula correct?  C? Patrick Is this formula correct?
353  cph Yes, but I rewrote it.  cph Yes, but I rewrote it.
354  cph Also, the KappaR? need the index and subscript k!  cph Also, the kappaR? need the index and subscript k!
355           kkey = (ikey-1)*Nr + k           kkey = (itdkey-1)*Nr + k
356  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
357    
358  C--       km1    Points to level above k (=k-1)  C--       km1    Points to level above k (=k-1)
# Line 521  C--       kDown  Cycles through 2,1 to p Line 368  C--       kDown  Cycles through 2,1 to p
368            jMin = 1-OLy            jMin = 1-OLy
369            jMax = sNy+OLy            jMax = sNy+OLy
370    
371              kp1Msk=1.
372              IF (k.EQ.Nr) kp1Msk=0.
373               DO j=1-Oly,sNy+Oly
374                DO i=1-Olx,sNx+Olx
375                 rTransKp1(i,j) = kp1Msk*rTrans(i,j)
376                ENDDO
377               ENDDO
378    #ifdef ALLOW_AUTODIFF_TAMC
379    CADJ STORE rTransKp1(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
380    #endif
381    
382  C--       Get temporary terms used by tendency routines  C--       Get temporary terms used by tendency routines
383            CALL CALC_COMMON_FACTORS (            CALL CALC_COMMON_FACTORS (
384       I         bi,bj,iMin,iMax,jMin,jMax,k,       I         bi,bj,iMin,iMax,jMin,jMax,k,
385       O         xA,yA,uTrans,vTrans,rTrans,maskUp,       O         xA,yA,uTrans,vTrans,rTrans,maskUp,
386       I         myThid)       I         myThid)
387    
388              IF (k.EQ.1) THEN
389    C- Surface interface :
390               DO j=1-Oly,sNy+Oly
391                DO i=1-Olx,sNx+Olx
392                 rTrans(i,j) = 0.
393                ENDDO
394               ENDDO
395              ELSE
396    C- Interior interface :
397               DO j=1-Oly,sNy+Oly
398                DO i=1-Olx,sNx+Olx
399                 rTrans(i,j) = rTrans(i,j)*maskC(i,j,k-1,bi,bj)
400                ENDDO
401               ENDDO
402              ENDIF
403    
404  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
405    
406  C--   Residual transp = Bolus transp + Eulerian transp  C--   Residual transp = Bolus transp + Eulerian transp
407            IF (useGMRedi) THEN            IF (useGMRedi) THEN
408              CALL GMREDI_CALC_UVFLOW(              CALL GMREDI_CALC_UVFLOW(
# Line 535  C--   Residual transp = Bolus transp + E Line 410  C--   Residual transp = Bolus transp + E
410              IF (K.GE.2) CALL GMREDI_CALC_WFLOW(              IF (K.GE.2) CALL GMREDI_CALC_WFLOW(
411       &                    rTrans, bi, bj, k, myThid)       &                    rTrans, bi, bj, k, myThid)
412            ENDIF            ENDIF
 #endif /* ALLOW_GMREDI */  
413    
414  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
415  CADJ STORE KappaRT(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE rTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
416  CADJ STORE KappaRS(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte  #ifdef GM_BOLUS_ADVEC
417    CADJ STORE uTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
418    CADJ STORE vTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
419    #endif
420  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
421    
422    #endif /* ALLOW_GMREDI */
423    
424  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL
425  C--      Calculate the total vertical diffusivity  C--      Calculate the total vertical diffusivity
426           CALL CALC_DIFFUSIVITY(            IF ( .NOT.implicitDiffusion ) THEN
427       I        bi,bj,iMin,iMax,jMin,jMax,k,              CALL CALC_DIFFUSIVITY(
428       I        maskUp,       I          bi,bj,iMin,iMax,jMin,jMax,k,
429       O        KappaRT,KappaRS,       I          maskUp,
430       I        myThid)       O          kappaRT,kappaRS,
431         I          myThid)
432              ENDIF
433    # ifdef ALLOW_AUTODIFF_TAMC
434    CADJ STORE kappaRT(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
435    CADJ STORE kappaRS(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
436    # endif /* ALLOW_AUTODIFF_TAMC */
437  #endif  #endif
438    
439            iMin = 1-OLx+2            iMin = 1-OLx+2
# Line 557  C--      Calculate the total vertical di Line 442  C--      Calculate the total vertical di
442            jMax = sNy+OLy-1            jMax = sNy+OLy-1
443    
444  C--      Calculate active tracer tendencies (gT,gS,...)  C--      Calculate active tracer tendencies (gT,gS,...)
445  C        and step forward storing result in gTnm1, gSnm1, etc.  C        and step forward storing result in gT, gS, etc.
446    #ifndef ALLOW_OFFLINE
447           IF ( tempStepping ) THEN           IF ( tempStepping ) THEN
448             CALL CALC_GT(             CALL CALC_GT(
449       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
450       I         xA,yA,uTrans,vTrans,rTrans,maskUp,       I         xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,
451       I         KappaRT,       I         kappaRT,
452       U         fVerT,       U         fVerT,
453       I         myTime,myIter,myThid)       I         myTime,myIter,myThid)
454             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
# Line 570  C        and step forward storing result Line 456  C        and step forward storing result
456       I         theta, gT,       I         theta, gT,
457       I         myIter, myThid)       I         myIter, myThid)
458           ENDIF           ENDIF
459    #endif
460    
461    #ifndef ALLOW_OFFLINE
462           IF ( saltStepping ) THEN           IF ( saltStepping ) THEN
463             CALL CALC_GS(             CALL CALC_GS(
464       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
465       I         xA,yA,uTrans,vTrans,rTrans,maskUp,       I         xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,
466       I         KappaRS,       I         kappaRS,
467       U         fVerS,       U         fVerS,
468       I         myTime,myIter,myThid)       I         myTime,myIter,myThid)
469             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
# Line 582  C        and step forward storing result Line 471  C        and step forward storing result
471       I         salt, gS,       I         salt, gS,
472       I         myIter, myThid)       I         myIter, myThid)
473           ENDIF           ENDIF
 #ifdef ALLOW_PASSIVE_TRACER  
          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)  
            CALL TIMESTEP_TRACER(  
      I         bi,bj,iMin,iMax,jMin,jMax,k,tracerAdvScheme,  
      I         Tr1, gTr1,  
      I         myIter,myThid)  
          ENDIF  
474  #endif  #endif
475  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
476           IF ( usePTRACERS ) THEN           IF ( usePTRACERS ) THEN
477             CALL PTRACERS_INTEGERATE(             IF ( .NOT.implicitDiffusion ) THEN
478                 CALL PTRACERS_CALC_DIFF(
479         I            bi,bj,iMin,iMax,jMin,jMax,k,
480         I            maskUp,
481         O            kappaRTr,
482         I            myThid)
483               ENDIF
484    # ifdef ALLOW_AUTODIFF_TAMC
485    CADJ STORE kappaRTr(:,:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
486    # endif /* ALLOW_AUTODIFF_TAMC */
487               CALL PTRACERS_INTEGRATE(
488       I         bi,bj,k,       I         bi,bj,k,
489       I         xA,yA,uTrans,vTrans,rTrans,maskUp,       I         xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,
490       X         KappaRS,       U         fVerP,
491         I         kappaRTr,
492       I         myIter,myTime,myThid)       I         myIter,myTime,myThid)
493           ENDIF           ENDIF
494  #endif /* ALLOW_PTRACERS */  #endif /* ALLOW_PTRACERS */
# Line 614  C--      Apply open boundary conditions Line 501  C--      Apply open boundary conditions
501  #endif   /* ALLOW_OBCS */  #endif   /* ALLOW_OBCS */
502    
503  C--      Freeze water  C--      Freeze water
504           IF (allowFreezing) THEN  C  this bit of code is left here for backward compatibility.
505    C  freezing at surface level has been moved to FORWARD_STEP
506    #ifndef ALLOW_OFFLINE
507             IF ( useOldFreezing .AND. .NOT. useSEAICE
508         &       .AND. .NOT.(useThSIce.AND.k.EQ.1) ) THEN
509  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
510  CADJ STORE gT(:,:,k,bi,bj) = comlev1_bibj_k  CADJ STORE gT(:,:,k,bi,bj) = comlev1_bibj_k
511  CADJ &   , key = kkey, byte = isbyte  CADJ &   , key = kkey, byte = isbyte
512  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
513              CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )              CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )
514           END IF           ENDIF
515    #endif
516    
517  C--     end of thermodynamic k loop (Nr:1)  C--     end of thermodynamic k loop (Nr:1)
518          ENDDO          ENDDO
519    
520    C       All explicit advection/diffusion/sources should now be
521  #ifdef ALLOW_AUTODIFF_TAMC  C       done. The updated tracer field is in gPtr. Accumalate
522  C? Patrick? What about this one?  C       explicit tendency and also reset gPtr to initial tracer
523  cph Keys iikey and idkey dont seem to be needed  C       field for implicit matrix calculation
524  cph since storing occurs on different tape for each  
525  cph impldiff call anyways.  #ifdef ALLOW_MATRIX
526  cph Thus, common block comlev1_impl isnt needed either.          IF (useMATRIX)
527  cph Storing below needed in the case useGMREDI.       &    CALL MATRIX_STORE_TENDENCY_EXP( bi, bj, myTime, myIter, myThid )
528          iikey = (ikey-1)*maximpl  #endif
529  #endif /* ALLOW_AUTODIFF_TAMC */  
530            iMin = 1
531  C--     Implicit diffusion          iMax = sNx
532          IF (implicitDiffusion) THEN          jMin = 1
533            jMax = sNy
534           IF (tempStepping) THEN  
535    C--     Implicit vertical advection & diffusion
536    #ifndef ALLOW_OFFLINE
537            IF ( tempStepping .AND. implicitDiffusion ) THEN
538              CALL CALC_3D_DIFFUSIVITY(
539         I         bi,bj,iMin,iMax,jMin,jMax,
540         I         GAD_TEMPERATURE, useGMredi, useKPP,
541         O         kappaRk,
542         I         myThid)
543            ENDIF
544    #ifdef INCLUDE_IMPLVERTADV_CODE
545            IF ( tempImplVertAdv ) THEN
546              CALL GAD_IMPLICIT_R(
547         I         tempImplVertAdv, tempAdvScheme, GAD_TEMPERATURE,
548         I         kappaRk, wVel, theta,
549         U         gT,
550         I         bi, bj, myTime, myIter, myThid )
551            ELSEIF ( tempStepping .AND. implicitDiffusion ) THEN
552    #else /* INCLUDE_IMPLVERTADV_CODE */
553            IF     ( tempStepping .AND. implicitDiffusion ) THEN
554    #endif /* INCLUDE_IMPLVERTADV_CODE */
555  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
556              idkey = iikey + 1  CADJ STORE kappaRk(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
557  CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
558  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
559              CALL IMPLDIFF(            CALL IMPLDIFF(
560       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
561       I         deltaTtracer, KappaRT, recip_HFacC,       I         GAD_TEMPERATURE, kappaRk, recip_hFacC,
562       U         gT,       U         gT,
563       I         myThid )       I         myThid )
564            ENDIF
565    #endif /* ndef ALLOW_OFFLINE */
566    
567    #ifdef ALLOW_TIMEAVE
568    #ifndef MINIMAL_TAVE_OUTPUT
569            useVariableK = useKPP .OR. usePP81 .OR. useMY82 .OR. useGGL90
570         &       .OR. useGMredi .OR. ivdc_kappa.NE.0.
571            IF (taveFreq.GT.0. .AND. useVariableK ) THEN
572             IF (implicitDiffusion) THEN
573               CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, gT, kappaRk,
574         I                        Nr, 3, deltaTclock, bi, bj, myThid)
575    c        ELSE
576    c          CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, theta, kappaRT,
577    c    I                        Nr, 3, deltaTclock, bi, bj, myThid)
578           ENDIF           ENDIF
579            ENDIF
580    #endif /* ndef MINIMAL_TAVE_OUTPUT */
581    #endif /* ALLOW_TIMEAVE */
582    
583           IF (saltStepping) THEN  #ifndef ALLOW_OFFLINE
584            IF ( saltStepping .AND. implicitDiffusion ) THEN
585              CALL CALC_3D_DIFFUSIVITY(
586         I         bi,bj,iMin,iMax,jMin,jMax,
587         I         GAD_SALINITY, useGMredi, useKPP,
588         O         kappaRk,
589         I         myThid)
590            ENDIF
591    
592    #ifdef INCLUDE_IMPLVERTADV_CODE
593            IF ( saltImplVertAdv ) THEN
594              CALL GAD_IMPLICIT_R(
595         I         saltImplVertAdv, saltAdvScheme, GAD_SALINITY,
596         I         kappaRk, wVel, salt,
597         U         gS,
598         I         bi, bj, myTime, myIter, myThid )
599            ELSEIF ( saltStepping .AND. implicitDiffusion ) THEN
600    #else /* INCLUDE_IMPLVERTADV_CODE */
601            IF     ( saltStepping .AND. implicitDiffusion ) THEN
602    #endif /* INCLUDE_IMPLVERTADV_CODE */
603  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
604           idkey = iikey + 2  CADJ STORE kappaRk(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
605  CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
606  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
607              CALL IMPLDIFF(            CALL IMPLDIFF(
608       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
609       I         deltaTtracer, KappaRS, recip_HFacC,       I         GAD_SALINITY, kappaRk, recip_hFacC,
610       U         gS,       U         gS,
611       I         myThid )       I         myThid )
612           ENDIF          ENDIF
   
 #ifdef ALLOW_PASSIVE_TRACER  
          IF (tr1Stepping) THEN  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE gTr1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, 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  
613  #endif  #endif
614    
615  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
616  C Vertical diffusion (implicit) for passive tracers          IF     ( usePTRACERS ) THEN
617           IF ( usePTRACERS ) THEN  C--     Vertical advection/diffusion (implicit) for passive tracers
618             CALL PTRACERS_IMPLDIFF( bi,bj,KappaRS,myThid )             CALL PTRACERS_IMPLICIT(
619           ENDIF       U                             kappaRk,
620         I                             bi, bj, myTime, myIter, myThid )
621            ENDIF
622  #endif /* ALLOW_PTRACERS */  #endif /* ALLOW_PTRACERS */
623    
624  #ifdef   ALLOW_OBCS  #ifdef   ALLOW_OBCS
625  C--      Apply open boundary conditions  C--      Apply open boundary conditions
626           IF (useOBCS) THEN          IF ( ( implicitDiffusion
627         &    .OR. tempImplVertAdv
628         &    .OR. saltImplVertAdv
629         &       ) .AND. useOBCS     ) THEN
630             DO K=1,Nr             DO K=1,Nr
631               CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )               CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )
632             ENDDO             ENDDO
633           END IF          ENDIF
634  #endif   /* ALLOW_OBCS */  #endif   /* ALLOW_OBCS */
635    
636  C--     End If implicitDiffusion  #ifdef ALLOW_TIMEAVE
637            IF ( taveFreq.GT. 0. _d 0 .AND. fluidIsWater ) THEN
638              CALL TIMEAVE_SURF_FLUX( bi, bj, myTime, myIter, myThid)
639          ENDIF          ENDIF
640    #ifndef MINIMAL_TAVE_OUTPUT
641            IF (taveFreq.GT.0. .AND. ivdc_kappa.NE.0.) THEN
642              CALL TIMEAVE_CUMULATE(ConvectCountTave, IVDConvCount,
643         I                           Nr, deltaTclock, bi, bj, myThid)
644            ENDIF
645    #endif /* ndef MINIMAL_TAVE_OUTPUT */
646    #endif /* ALLOW_TIMEAVE */
647    
648    #ifdef ALLOW_DIAGNOSTICS
649            IF ( usediagnostics ) CALL DIAGNOSTICS_FILL_SURF_FLUX(myThid)
650    #endif
651    
652  #endif /* SINGLE_LAYER_MODE */  #endif /* SINGLE_LAYER_MODE */
653    
654  Ccs-  C--   end bi,bj loops.
655         ENDDO         ENDDO
656        ENDDO        ENDDO
657    
658  #ifdef ALLOW_AIM  #ifdef ALLOW_DEBUG
       IF ( useAIM ) THEN  
        CALL AIM_AIM2DYN_EXCHANGES( myTime, myThid )  
       ENDIF  
 #endif /* ALLOW_AIM */  
       IF ( staggerTimeStep ) THEN  
        IF ( useAIM .OR. useCubedSphereExchange ) THEN  
          IF (tempStepping) _EXCH_XYZ_R8(gT,myThid)  
          IF (saltStepping) _EXCH_XYZ_R8(gS,myThid)  
        ELSEIF ( useGMRedi .AND. Oly.LT.4 ) THEN  
 c        .AND. GM_AdvForm .AND. .NOT.GM_AdvSeparate ) THEN  
          IF (tempMultiDimAdvec) _EXCH_XYZ_R8(gT,myThid)  
          IF (saltMultiDimAdvec) _EXCH_XYZ_R8(gS,myThid)  
        ENDIF    
       ENDIF    
   
 #ifndef DISABLE_DEBUGMODE  
659        If (debugMode) THEN        If (debugMode) THEN
660         CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (THERMODYNAMICS)',myThid)
661         CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (THERMODYNAMICS)',myThid)
662         CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (THERMODYNAMICS)',myThid)
663         CALL DEBUG_STATS_RL(Nr,theta,'Theta (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,theta,'Theta (THERMODYNAMICS)',myThid)
664         CALL DEBUG_STATS_RL(Nr,salt,'Salt (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,salt,'Salt (THERMODYNAMICS)',myThid)
665         CALL DEBUG_STATS_RL(Nr,Gt,'Gt (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,gT,'Gt (THERMODYNAMICS)',myThid)
666         CALL DEBUG_STATS_RL(Nr,Gs,'Gs (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,gS,'Gs (THERMODYNAMICS)',myThid)
667         CALL DEBUG_STATS_RL(Nr,GtNm1,'GtNm1 (THERMODYNAMICS)',myThid)  #ifndef ALLOW_ADAMSBASHFORTH_3
668         CALL DEBUG_STATS_RL(Nr,GsNm1,'GsNm1 (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,gtNm1,'GtNm1 (THERMODYNAMICS)',myThid)
669           CALL DEBUG_STATS_RL(Nr,gsNm1,'GsNm1 (THERMODYNAMICS)',myThid)
670    #endif
671  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
672         IF ( usePTRACERS ) THEN         IF ( usePTRACERS ) THEN
673           CALL PTRACERS_DEBUG(myThid)           CALL PTRACERS_DEBUG(myThid)
674         ENDIF         ENDIF
675  #endif /* ALLOW_PTRACERS */  #endif /* ALLOW_PTRACERS */
676        ENDIF        ENDIF
677    #endif /* ALLOW_DEBUG */
678    
679    #ifdef ALLOW_DEBUG
680             IF ( debugLevel .GE. debLevB )
681         &    CALL DEBUG_LEAVE('THERMODYNAMICS',myThid)
682  #endif  #endif
683    
684    #endif /* ALLOW_GENERIC_ADVDIFF */
685    
686        RETURN        RETURN
687        END        END
Legend:
Removed from v.1.24  
changed lines
  Added in v.1.91
  ViewVC Help
Powered by ViewVC 1.1.22