/[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.112 by heimbach, Sat Jan 6 22:30:17 2007 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
8    # ifdef ALLOW_GMREDI
9    #  include "GMREDI_OPTIONS.h"
10    # endif
11    # ifdef ALLOW_KPP
12    #  include "KPP_OPTIONS.h"
13    # endif
14    #endif /* ALLOW_AUTODIFF_TAMC */
15    
16  CBOP  CBOP
17  C     !ROUTINE: THERMODYNAMICS  C     !ROUTINE: THERMODYNAMICS
18  C     !INTERFACE:  C     !INTERFACE:
19        SUBROUTINE THERMODYNAMICS(myTime, myIter, myThid)        SUBROUTINE THERMODYNAMICS(myTime, myIter, myThid)
20  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
21  C     *==========================================================*  C     *==========================================================*
22  C     | SUBROUTINE THERMODYNAMICS                                  C     | SUBROUTINE THERMODYNAMICS
23  C     | o Controlling routine for the prognostic part of the        C     | o Controlling routine for the prognostic part of the
24  C     |   thermo-dynamics.                                          C     |   thermo-dynamics.
25  C     *===========================================================  C     *===========================================================
26  C     | The algorithm...  C     | The algorithm...
27  C     |  C     |
# Line 67  C     == Global variables === Line 77  C     == Global variables ===
77  #include "PARAMS.h"  #include "PARAMS.h"
78  #include "DYNVARS.h"  #include "DYNVARS.h"
79  #include "GRID.h"  #include "GRID.h"
80    #ifdef ALLOW_GENERIC_ADVDIFF
81  #include "GAD.h"  #include "GAD.h"
 #ifdef ALLOW_PASSIVE_TRACER  
 #include "TR1.h"  
82  #endif  #endif
83    #ifdef ALLOW_PTRACERS
84    #include "PTRACERS_SIZE.h"
85    #include "PTRACERS.h"
86    #endif
87    #ifdef ALLOW_TIMEAVE
88    #include "TIMEAVE_STATV.h"
89    #endif
90    
91  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
92  # include "tamc.h"  # include "tamc.h"
93  # include "tamc_keys.h"  # include "tamc_keys.h"
94  # include "FFIELDS.h"  # include "FFIELDS.h"
95    # include "EOS.h"
96  # ifdef ALLOW_KPP  # ifdef ALLOW_KPP
97  #  include "KPP.h"  #  include "KPP.h"
98  # endif  # endif
99  # ifdef ALLOW_GMREDI  # ifdef ALLOW_GMREDI
100  #  include "GMREDI.h"  #  include "GMREDI.h"
101  # endif  # endif
102    # ifdef ALLOW_EBM
103    #  include "EBM.h"
104    # endif
105  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
 #ifdef ALLOW_TIMEAVE  
 #include "TIMEAVE_STATV.h"  
 #endif  
106    
107  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
108  C     == Routine arguments ==  C     == Routine arguments ==
# Line 95  C     myThid - Thread number for this in Line 113  C     myThid - Thread number for this in
113        INTEGER myIter        INTEGER myIter
114        INTEGER myThid        INTEGER myThid
115    
116    #ifdef ALLOW_GENERIC_ADVDIFF
117  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
118  C     == Local variables  C     == Local variables
119  C     xA, yA                 - Per block temporaries holding face areas  C     xA, yA                 - Per block temporaries holding face areas
120  C     uTrans, vTrans, rTrans - Per block temporaries holding flow  C     uFld, vFld, wFld       - Local copy of velocity field (3 components)
121  C                              transport  C     uTrans, vTrans, rTrans - Per block temporaries holding flow transport
122  C                              o uTrans: Zonal transport  C                              o uTrans: Zonal transport
123  C                              o vTrans: Meridional transport  C                              o vTrans: Meridional transport
124  C                              o rTrans: Vertical transport  C                              o rTrans: Vertical transport
125    C     rTransKp1                o vertical volume transp. at interface k+1
126  C     maskUp                   o maskUp: land/water mask for W points  C     maskUp                   o maskUp: land/water mask for W points
127  C     fVer[STUV]               o fVer: Vertical flux term - note fVer  C     fVer[STUV]               o fVer: Vertical flux term - note fVer
128  C                                      is "pipelined" in the vertical  C                                      is "pipelined" in the vertical
129  C                                      so we need an fVer for each  C                                      so we need an fVer for each
130  C                                      variable.  C                                      variable.
131  C     rhoK, rhoKM1   - Density at current level, and level above  C     kappaRT,       - Total diffusion in vertical at level k, for T and S
132  C     phiHyd         - Hydrostatic part of the potential phiHydi.  C     kappaRS          (background + spatially varying, isopycnal term).
133  C                      In z coords phiHydiHyd is the hydrostatic  C     kappaRTr       - Total diffusion in vertical at level k,
134  C                      Potential (=pressure/rho0) anomaly  C                      for each passive Tracer
135  C                      In p coords phiHydiHyd is the geopotential  C     kappaRk        - Total diffusion in vertical, all levels, 1 tracer
136  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).  
137  C     iMin, iMax     - Ranges and sub-block indices on which calculations  C     iMin, iMax     - Ranges and sub-block indices on which calculations
138  C     jMin, jMax       are applied.  C     jMin, jMax       are applied.
139  C     bi, bj  C     bi, bj
# Line 126  C     kDown, km1       are switched with Line 142  C     kDown, km1       are switched with
142  C                      index into fVerTerm.  C                      index into fVerTerm.
143        _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
144        _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
145          _RL uFld    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
146          _RL vFld    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
147          _RL wFld    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
148        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
149        _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
150        _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
151          _RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
152        _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
153        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
154        _RL fVerS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
155        _RL fVerTr1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL kappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
156        _RL phiHyd  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL kappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
157        _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  #ifdef ALLOW_PTRACERS
158        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fVerP   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2,PTRACERS_num)
159        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL kappaRTr(1-Olx:sNx+Olx,1-Oly:sNy+Oly,PTRACERS_num)
160        _RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  #endif
161        _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL kappaRk (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
       _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)  
       _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
       _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
       _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
 C     This is currently used by IVDC and Diagnostics  
       _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
162        INTEGER iMin, iMax        INTEGER iMin, iMax
163        INTEGER jMin, jMax        INTEGER jMin, jMax
164        INTEGER bi, bj        INTEGER bi, bj
165        INTEGER i, j        INTEGER i, j
166        INTEGER k, km1, kup, kDown        INTEGER k, km1, kup, kDown
167    #ifdef ALLOW_ADAMSBASHFORTH_3
168          INTEGER iterNb, m1, m2
169    #endif
170    #ifdef ALLOW_TIMEAVE
171          LOGICAL useVariableK
172    #endif
173    #ifdef ALLOW_PTRACERS
174          INTEGER iTracer, ip
175    #endif
176    
177  CEOP  CEOP
178    
179    #ifdef ALLOW_DEBUG
180             IF ( debugLevel .GE. debLevB )
181         &    CALL DEBUG_ENTER('THERMODYNAMICS',myThid)
182    #endif
183    
184  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
185  C--   dummy statement to end declaration part  C--   dummy statement to end declaration part
186        ikey = 1        ikey = 1
187          itdkey = 1
188  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
189    
 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  
   
   
190  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
191  C--   HPF directive to help TAMC  C--   HPF directive to help TAMC
192  CHPF$ INDEPENDENT  CHPF$ INDEPENDENT
193  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
194    
195    C-- Compute correction at the surface for Lin Free Surf.
196          IF (linFSConserveTr)
197         &   CALL CALC_WSURF_TR(theta,salt,wVel,
198         &                      myIter,myTime,myThid)
199    
200        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
201    
202  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
203  C--    HPF directive to help TAMC  C--    HPF directive to help TAMC
204  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS
205  CHPF$&                  ,phiHyd,utrans,vtrans,xA,yA  CHPF$&                  ,utrans,vtrans,xA,yA
206  CHPF$&                  ,KappaRT,KappaRS  CHPF$&                  ,kappaRT,kappaRS
207  CHPF$&                  )  CHPF$&                  )
208    # ifdef ALLOW_PTRACERS
209    CHPF$  INDEPENDENT, NEW (fVerP,kappaRTr)
210    # endif
211  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
212    
213         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
# Line 208  CHPF$&                  ) Line 220  CHPF$&                  )
220            act3 = myThid - 1            act3 = myThid - 1
221            max3 = nTx*nTy            max3 = nTx*nTy
222            act4 = ikey_dynamics - 1            act4 = ikey_dynamics - 1
223            ikey = (act1 + 1) + act2*max1            itdkey = (act1 + 1) + act2*max1
224       &                      + act3*max1*max2       &                      + act3*max1*max2
225       &                      + act4*max1*max2*max3       &                      + act4*max1*max2*max3
226  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
227    
228  C--     Set up work arrays that need valid initial values  C--   Set up work arrays with valid (i.e. not NaN) values
229    C     These inital values do not alter the numerical results. They
230    C     just ensure that all memory references are to valid floating
231    C     point numbers. This prevents spurious hardware signals due to
232    C     uninitialised but inert locations.
233    
234          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
235           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
236              xA(i,j)        = 0. _d 0
237              yA(i,j)        = 0. _d 0
238              uTrans(i,j)    = 0. _d 0
239              vTrans(i,j)    = 0. _d 0
240            rTrans (i,j)   = 0. _d 0            rTrans (i,j)   = 0. _d 0
241              rTransKp1(i,j) = 0. _d 0
242            fVerT  (i,j,1) = 0. _d 0            fVerT  (i,j,1) = 0. _d 0
243            fVerT  (i,j,2) = 0. _d 0            fVerT  (i,j,2) = 0. _d 0
244            fVerS  (i,j,1) = 0. _d 0            fVerS  (i,j,1) = 0. _d 0
245            fVerS  (i,j,2) = 0. _d 0            fVerS  (i,j,2) = 0. _d 0
246            fVerTr1(i,j,1) = 0. _d 0            kappaRT(i,j)   = 0. _d 0
247            fVerTr1(i,j,2) = 0. _d 0            kappaRS(i,j)   = 0. _d 0
248           ENDDO           ENDDO
249          ENDDO          ENDDO
250    
# Line 230  C--     Set up work arrays that need val Line 252  C--     Set up work arrays that need val
252           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
253            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
254  C This is currently also used by IVDC and Diagnostics  C This is currently also used by IVDC and Diagnostics
255             ConvectCount(i,j,k) = 0.             kappaRk(i,j,k)    = 0. _d 0
256             KappaRT(i,j,k) = 0. _d 0  C- tracer tendency needs to be set to zero (moved here from gad_calc_rhs):
257             KappaRS(i,j,k) = 0. _d 0             gT(i,j,k,bi,bj)   = 0. _d 0
258  #ifdef ALLOW_AUTODIFF_TAMC             gS(i,j,k,bi,bj)   = 0. _d 0
            gT(i,j,k,bi,bj) = 0. _d 0  
            gS(i,j,k,bi,bj) = 0. _d 0  
 #ifdef ALLOW_PASSIVE_TRACER  
            gTr1(i,j,k,bi,bj) = 0. _d 0  
 #endif  
 #endif  
259            ENDDO            ENDDO
260           ENDDO           ENDDO
261          ENDDO          ENDDO
262    
263          iMin = 1-OLx+1  #ifdef ALLOW_PTRACERS
264          iMax = sNx+OLx          IF ( usePTRACERS ) THEN
265          jMin = 1-OLy+1           DO ip=1,PTRACERS_num
266          jMax = sNy+OLy             DO j=1-OLy,sNy+OLy
267                DO i=1-OLx,sNx+OLx
268                 fVerP  (i,j,1,ip) = 0. _d 0
269  #ifdef ALLOW_AUTODIFF_TAMC               fVerP  (i,j,2,ip) = 0. _d 0
270  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte               kappaRTr(i,j,ip)  = 0. _d 0
271  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte              ENDDO
272  #ifdef ALLOW_KPP             ENDDO
273  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte           ENDDO
274  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  C-      set tracer tendency to zero:
275  #endif           DO iTracer=1,PTRACERS_num
 #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--       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  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 cph avoids recomputation of integrate_for_w  
 CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, 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,  
      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  
276            DO k=1,Nr            DO k=1,Nr
277              CALL GMREDI_CALC_TENSOR(             DO j=1-OLy,sNy+OLy
278       I             bi, bj, iMin, iMax, jMin, jMax, k,              DO i=1-OLx,sNx+OLx
279       I             sigmaX, sigmaY, sigmaR,               gPTr(i,j,k,bi,bj,itracer) = 0. _d 0
280       I             myThid )              ENDDO
281            ENDDO             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 )  
282            ENDDO            ENDDO
283  #endif /* ALLOW_AUTODIFF_TAMC */           ENDDO
284          ENDIF          ENDIF
285    #endif
286    
287  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_ADAMSBASHFORTH_3
288  CADJ STORE Kwx(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte  C-      Apply AB on T,S :
289  CADJ STORE Kwy(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte          iterNb = myIter
290  CADJ STORE Kwz(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte          IF (staggerTimeStep) iterNb = myIter - 1
291  #endif /* ALLOW_AUTODIFF_TAMC */          m1 = 1 + MOD(iterNb+1,2)
292            m2 = 1 + MOD( iterNb ,2)
293  #endif  /* ALLOW_GMREDI */  C       compute T^n+1/2 (stored in gtNm) extrapolating T forward in time
294            IF ( AdamsBashforth_T ) CALL ADAMS_BASHFORTH3(
295  #ifdef  ALLOW_KPP       I                                  bi, bj, 0,
296  C--     Compute KPP mixing coefficients       U                                  theta, gtNm,
297          IF (useKPP) THEN       I                                  tempStartAB, iterNb, myThid )
298            CALL KPP_CALC(  C       compute S^n+1/2 (stored in gsNm) extrapolating S forward in time
299       I                  bi, bj, myTime, myThid )          IF ( AdamsBashforth_S ) CALL ADAMS_BASHFORTH3(
300  #ifdef ALLOW_AUTODIFF_TAMC       I                                  bi, bj, 0,
301          ELSE       U                                  salt, gsNm,
302            CALL KPP_CALC_DUMMY(       I                                  saltStartAB, iterNb, myThid )
303       I                  bi, bj, myTime, myThid )  #endif /* ALLOW_ADAMSBASHFORTH_3 */
304  #endif /* ALLOW_AUTODIFF_TAMC */  
305          ENDIF  c       iMin = 1-OLx
306    c       iMax = sNx+OLx
307    c       jMin = 1-OLy
308    c       jMax = sNy+OLy
309    
310  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
311  CADJ STORE KPPghat   (:,:,:,bi,bj)  cph avoids recomputation of integrate_for_w
312  CADJ &   , KPPdiffKzT(:,:,:,bi,bj)  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
 CADJ &   , KPPdiffKzS(:,:,:,bi,bj)  
 CADJ &   , KPPfrac   (:,:  ,bi,bj)  
 CADJ &                 = comlev1_bibj, key=ikey, byte=isbyte  
313  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
314    
315  #endif  /* ALLOW_KPP */  C--     Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h
316    C--     MOST of THERMODYNAMICS will be disabled
317    #ifndef SINGLE_LAYER_MODE
318    
319  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
320  CADJ STORE KappaRT(:,:,:)     = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
321  CADJ STORE KappaRS(:,:,:)     = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
322  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
323  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
324  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  # ifdef ALLOW_DEPTH_CONTROL
325  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE gtnm1(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
326  #ifdef ALLOW_PASSIVE_TRACER  CADJ STORE gsnm1(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
327  CADJ STORE tr1  (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  # endif
 #endif  
328  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
329    
 #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 */  
   
330  #ifndef DISABLE_MULTIDIM_ADVECTION  #ifndef DISABLE_MULTIDIM_ADVECTION
331  C--     Some advection schemes are better calculated using a multi-dimensional  C--     Some advection schemes are better calculated using a multi-dimensional
332  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 338  C to be able to exclude this scheme to a
338  C recomputation. It *is* differentiable, if you need it.  C recomputation. It *is* differentiable, if you need it.
339  C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to  C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to
340  C disable this section of code.  C disable this section of code.
341          IF (multiDimAdvection) THEN          IF (tempMultiDimAdvec) THEN
342           IF (tempStepping .AND.  #ifdef ALLOW_DEBUG
343       &       tempAdvScheme.NE.ENUM_CENTERED_2ND .AND.            IF ( debugLevel .GE. debLevB )
344       &       tempAdvScheme.NE.ENUM_UPWIND_3RD .AND.       &     CALL DEBUG_CALL('GAD_ADVECTION',myThid)
345       &       tempAdvScheme.NE.ENUM_CENTERED_4TH ) THEN  #endif
346            CALL GAD_ADVECTION(bi,bj,tempAdvScheme,GAD_TEMPERATURE,            CALL GAD_ADVECTION(
347       U                      theta,gT,       I             tempImplVertAdv, tempAdvScheme, tempVertAdvScheme,
348       I                      myTime,myIter,myThid)       I             GAD_TEMPERATURE,
349           ENDIF       I             uVel, vVel, wVel, theta,
350           IF (saltStepping .AND.       O             gT,
351       &       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  
352          ENDIF          ENDIF
353            IF (saltMultiDimAdvec) THEN
354    #ifdef ALLOW_DEBUG
355              IF ( debugLevel .GE. debLevB )
356         &     CALL DEBUG_CALL('GAD_ADVECTION',myThid)
357    #endif
358              CALL GAD_ADVECTION(
359         I             saltImplVertAdv, saltAdvScheme, saltVertAdvScheme,
360         I             GAD_SALINITY,
361         I             uVel, vVel, wVel, salt,
362         O             gS,
363         I             bi,bj,myTime,myIter,myThid)
364            ENDIF
365    
366    C Since passive tracers are configurable separately from T,S we
367    C call the multi-dimensional method for PTRACERS regardless
368    C of whether multiDimAdvection is set or not.
369    #ifdef ALLOW_PTRACERS
370            IF ( usePTRACERS ) THEN
371    #ifdef ALLOW_DEBUG
372              IF ( debugLevel .GE. debLevB )
373         &     CALL DEBUG_CALL('PTRACERS_ADVECTION',myThid)
374    #endif
375             CALL PTRACERS_ADVECTION( bi,bj,myIter,myTime,myThid )
376            ENDIF
377    #endif /* ALLOW_PTRACERS */
378  #endif /* DISABLE_MULTIDIM_ADVECTION */  #endif /* DISABLE_MULTIDIM_ADVECTION */
379    
380    #ifdef ALLOW_DEBUG
381            IF ( debugLevel .GE. debLevB )
382         &    CALL DEBUG_MSG('ENTERING DOWNWARD K LOOP',myThid)
383    #endif
384    
385  C--     Start of thermodynamics loop  C--     Start of thermodynamics loop
386          DO k=Nr,1,-1          DO k=Nr,1,-1
387  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
388  C? Patrick Is this formula correct?  C? Patrick Is this formula correct?
389  cph Yes, but I rewrote it.  cph Yes, but I rewrote it.
390  cph Also, the KappaR? need the index and subscript k!  cph Also, the kappaR? need the index and subscript k!
391           kkey = (ikey-1)*Nr + k           kkey = (itdkey-1)*Nr + k
392  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
393    
394  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 404  C--       kDown  Cycles through 2,1 to p
404            jMin = 1-OLy            jMin = 1-OLy
405            jMax = sNy+OLy            jMax = sNy+OLy
406    
407  C--       Get temporary terms used by tendency routines            IF (k.EQ.Nr) THEN
408               DO j=1-Oly,sNy+Oly
409                DO i=1-Olx,sNx+Olx
410                 rTransKp1(i,j) = 0. _d 0
411                ENDDO
412               ENDDO
413              ELSE
414               DO j=1-Oly,sNy+Oly
415                DO i=1-Olx,sNx+Olx
416                 rTransKp1(i,j) = rTrans(i,j)
417                ENDDO
418               ENDDO
419              ENDIF
420    #ifdef ALLOW_AUTODIFF_TAMC
421    CADJ STORE rTransKp1(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
422    #endif
423    
424    C--       Get temporary terms used by tendency routines :
425    C-        Calculate horizontal "volume transport" through tracer cell face
426    C         anelastic: uTrans,vTrans are scaled by rhoFacC (~ mass transport)
427            CALL CALC_COMMON_FACTORS (            CALL CALC_COMMON_FACTORS (
428       I         bi,bj,iMin,iMax,jMin,jMax,k,       I         uVel, vVel,
429       O         xA,yA,uTrans,vTrans,rTrans,maskUp,       O         uFld, vFld, uTrans, vTrans, xA, yA,
430       I         myThid)       I         k,bi,bj, myThid )
431    
432    C-        Calculate vertical "volume transport" through tracer cell face
433              IF (k.EQ.1) THEN
434    C-         Surface interface :
435               DO j=1-Oly,sNy+Oly
436                DO i=1-Olx,sNx+Olx
437                 wFld(i,j)   = 0. _d 0
438                 maskUp(i,j) = 0. _d 0
439                 rTrans(i,j) = 0. _d 0
440                ENDDO
441               ENDDO
442              ELSE
443    C-         Interior interface :
444    C          anelastic: rTrans is scaled by rhoFacF (~ mass transport)
445               DO j=1-Oly,sNy+Oly
446                DO i=1-Olx,sNx+Olx
447                 wFld(i,j)   = wVel(i,j,k,bi,bj)
448                 maskUp(i,j) = maskC(i,j,k-1,bi,bj)*maskC(i,j,k,bi,bj)
449                 rTrans(i,j) = wFld(i,j)*rA(i,j,bi,bj)*maskUp(i,j)
450         &                              *deepFac2F(k)*rhoFacF(k)
451                ENDDO
452               ENDDO
453              ENDIF
454    
455  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_GMREDI
456  CADJ STORE KappaRT(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte  C--   Residual transp = Bolus transp + Eulerian transp
457  CADJ STORE KappaRS(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte            IF (useGMRedi) THEN
458  #endif /* ALLOW_AUTODIFF_TAMC */              CALL GMREDI_CALC_UVFLOW(
459         U                  uFld, vFld, uTrans, vTrans,
460         I                  k, bi, bj, myThid )
461                IF (K.GE.2) THEN
462                  CALL GMREDI_CALC_WFLOW(
463         U                  wFld, rTrans,
464         I                  k, bi, bj, myThid )
465                ENDIF
466              ENDIF
467    # ifdef ALLOW_AUTODIFF_TAMC
468    CADJ STORE rTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
469    CADJ STORE wfld(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
470    # ifdef GM_BOLUS_ADVEC
471    CADJ STORE ufld(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
472    CADJ STORE vfld(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
473    CADJ STORE uTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
474    CADJ STORE vTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
475    # endif
476    # endif /* ALLOW_AUTODIFF_TAMC */
477    #endif /* ALLOW_GMREDI */
478    
479  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL
480  C--      Calculate the total vertical diffusivity  C--      Calculate the total vertical diffusivity
481           CALL CALC_DIFFUSIVITY(            IF ( .NOT.implicitDiffusion ) THEN
482       I        bi,bj,iMin,iMax,jMin,jMax,k,              CALL CALC_DIFFUSIVITY(
483       I        maskUp,       I          bi,bj,iMin,iMax,jMin,jMax,k,
484       O        KappaRT,KappaRS,       I          maskUp,
485       I        myThid)       O          kappaRT,kappaRS,
486         I          myThid)
487              ENDIF
488    # ifdef ALLOW_AUTODIFF_TAMC
489    CADJ STORE kappaRT(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
490    CADJ STORE kappaRS(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
491    # endif /* ALLOW_AUTODIFF_TAMC */
492  #endif  #endif
493    
494            iMin = 1-OLx+2            iMin = 1-OLx+2
# Line 526  C--      Calculate the total vertical di Line 497  C--      Calculate the total vertical di
497            jMax = sNy+OLy-1            jMax = sNy+OLy-1
498    
499  C--      Calculate active tracer tendencies (gT,gS,...)  C--      Calculate active tracer tendencies (gT,gS,...)
500  C        and step forward storing result in gTnm1, gSnm1, etc.  C        and step forward storing result in gT, gS, etc.
501    C--
502    # ifdef ALLOW_AUTODIFF_TAMC
503    #  if ((defined NONLIN_FRSURF) || (defined ALLOW_DEPTH_CONTROL)) && (defined ALLOW_GMREDI)
504    CADJ STORE kux(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
505    CADJ STORE kvy(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
506    #   ifdef GM_EXTRA_DIAGONAL
507    CADJ STORE kuz(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
508    CADJ STORE kvz(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
509    #   endif
510    #  endif
511    # endif /* ALLOW_AUTODIFF_TAMC */
512    C
513    #ifdef ALLOW_AUTODIFF_TAMC
514    # if (defined NONLIN_FRSURF) || (defined ALLOW_DEPTH_CONTROL)
515    cph-test
516    CADJ STORE uFld(:,:), vFld(:,:), wFld(:,:)
517    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
518    CADJ STORE uTrans(:,:), vTrans(:,:)
519    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
520    CADJ STORE xA(:,:), yA(:,:)
521    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
522    # endif
523    #endif /* ALLOW_AUTODIFF_TAMC */
524    C
525           IF ( tempStepping ) THEN           IF ( tempStepping ) THEN
526             CALL CALC_GT(  #ifdef ALLOW_AUTODIFF_TAMC
527    CADJ STORE gTnm1(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
528    # if (defined NONLIN_FRSURF) || (defined ALLOW_DEPTH_CONTROL)
529    CADJ STORE fvert(:,:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
530    # endif
531    #endif /* ALLOW_AUTODIFF_TAMC */
532              CALL CALC_GT(
533       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
534       I         xA,yA,uTrans,vTrans,rTrans,maskUp,       I         xA, yA, maskUp, uFld, vFld, wFld,
535       I         KappaRT,       I         uTrans, vTrans, rTrans, rTransKp1,
536         I         kappaRT,
537       U         fVerT,       U         fVerT,
538       I         myTime,myIter,myThid)       I         myTime,myIter,myThid)
539    #ifdef ALLOW_ADAMSBASHFORTH_3
540              IF ( AdamsBashforth_T ) THEN
541               CALL TIMESTEP_TRACER(
542         I         bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,
543         I         gtNm(1-Olx,1-Oly,1,1,1,m2),
544         U         gT,
545         I         myIter, myThid)
546              ELSE
547    #endif
548             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
549       I         bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,       I         bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,
550       I         theta, gT,       I         theta,
551         U         gT,
552       I         myIter, myThid)       I         myIter, myThid)
553    #ifdef ALLOW_ADAMSBASHFORTH_3
554              ENDIF
555    #endif
556           ENDIF           ENDIF
557    
558           IF ( saltStepping ) THEN           IF ( saltStepping ) THEN
559             CALL CALC_GS(  #ifdef ALLOW_AUTODIFF_TAMC
560    CADJ STORE gSnm1(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
561    # if (defined NONLIN_FRSURF) || (defined ALLOW_DEPTH_CONTROL)
562    CADJ STORE fvers(:,:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
563    # endif
564    #endif /* ALLOW_AUTODIFF_TAMC */
565              CALL CALC_GS(
566       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
567       I         xA,yA,uTrans,vTrans,rTrans,maskUp,       I         xA, yA, maskUp, uFld, vFld, wFld,
568       I         KappaRS,       I         uTrans, vTrans, rTrans, rTransKp1,
569         I         kappaRS,
570       U         fVerS,       U         fVerS,
571       I         myTime,myIter,myThid)       I         myTime,myIter,myThid)
572    #ifdef ALLOW_ADAMSBASHFORTH_3
573              IF ( AdamsBashforth_S ) THEN
574             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
575       I         bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,       I         bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,
576       I         salt, gS,       I         gsNm(1-Olx,1-Oly,1,1,1,m2),
577         U         gS,
578       I         myIter, myThid)       I         myIter, myThid)
579           ENDIF            ELSE
580  #ifdef ALLOW_PASSIVE_TRACER  #endif
          IF ( tr1Stepping ) THEN  
            CALL CALC_GTR1(  
      I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,  
      I         xA,yA,uTrans,vTrans,rTrans,maskUp,  
      I         KappaRT,  
      U         fVerTr1,  
      I         myTime,myIter,myThid)  
581             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
582       I         bi,bj,iMin,iMax,jMin,jMax,k,tracerAdvScheme,       I         bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,
583       I         Tr1, gTr1,       I         salt,
584       I         myIter,myThid)       U         gS,
585           ENDIF       I         myIter, myThid)
586    #ifdef ALLOW_ADAMSBASHFORTH_3
587              ENDIF
588  #endif  #endif
589             ENDIF
590    
591    #ifdef ALLOW_PTRACERS
592             IF ( usePTRACERS ) THEN
593               IF ( .NOT.implicitDiffusion ) THEN
594                 CALL PTRACERS_CALC_DIFF(
595         I            bi,bj,iMin,iMax,jMin,jMax,k,
596         I            maskUp,
597         O            kappaRTr,
598         I            myThid)
599               ENDIF
600    # ifdef ALLOW_AUTODIFF_TAMC
601    CADJ STORE kappaRTr(:,:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
602    # endif /* ALLOW_AUTODIFF_TAMC */
603               CALL PTRACERS_INTEGRATE(
604         I          bi,bj,k,
605         I          xA, yA, maskUp, uFld, vFld, wFld,
606         I          uTrans, vTrans, rTrans, rTransKp1,
607         I          kappaRTr,
608         U          fVerP,
609         I          myTime,myIter,myThid)
610             ENDIF
611    #endif /* ALLOW_PTRACERS */
612    
613  #ifdef   ALLOW_OBCS  #ifdef   ALLOW_OBCS
614  C--      Apply open boundary conditions  C--      Apply open boundary conditions
# Line 574  C--      Apply open boundary conditions Line 618  C--      Apply open boundary conditions
618  #endif   /* ALLOW_OBCS */  #endif   /* ALLOW_OBCS */
619    
620  C--      Freeze water  C--      Freeze water
621           IF (allowFreezing) THEN  C  this bit of code is left here for backward compatibility.
622    C  freezing at surface level has been moved to FORWARD_STEP
623             IF ( useOldFreezing .AND. .NOT. useSEAICE
624         &       .AND. .NOT.(useThSIce.AND.k.EQ.1) ) THEN
625  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
626  CADJ STORE gT(:,:,k,bi,bj) = comlev1_bibj_k  CADJ STORE gT(:,:,k,bi,bj) = comlev1_bibj_k
627  CADJ &   , key = kkey, byte = isbyte  CADJ &   , key = kkey, byte = isbyte
628  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
629              CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )              CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )
630           END IF           ENDIF
631    
632  C--     end of thermodynamic k loop (Nr:1)  C--     end of thermodynamic k loop (Nr:1)
633          ENDDO          ENDDO
634    
635    C       All explicit advection/diffusion/sources should now be
636  #ifdef ALLOW_AUTODIFF_TAMC  C       done. The updated tracer field is in gPtr. Accumalate
637  C? Patrick? What about this one?  C       explicit tendency and also reset gPtr to initial tracer
638  cph Keys iikey and idkey dont seem to be needed  C       field for implicit matrix calculation
639  cph since storing occurs on different tape for each  
640  cph impldiff call anyways.  #ifdef ALLOW_MATRIX
641  cph Thus, common block comlev1_impl isnt needed either.          IF (useMATRIX)
642  cph Storing below needed in the case useGMREDI.       &    CALL MATRIX_STORE_TENDENCY_EXP(bi,bj, myTime,myIter,myThid)
643          iikey = (ikey-1)*maximpl  #endif
644  #endif /* ALLOW_AUTODIFF_TAMC */  
645            iMin = 1
646  C--     Implicit diffusion          iMax = sNx
647          IF (implicitDiffusion) THEN          jMin = 1
648            jMax = sNy
649           IF (tempStepping) THEN  
650    C--     Implicit vertical advection & diffusion
651            IF ( tempStepping .AND. implicitDiffusion ) THEN
652              CALL CALC_3D_DIFFUSIVITY(
653         I         bi,bj,iMin,iMax,jMin,jMax,
654         I         GAD_TEMPERATURE, useGMredi, useKPP,
655         O         kappaRk,
656         I         myThid)
657            ENDIF
658    #ifdef INCLUDE_IMPLVERTADV_CODE
659            IF ( tempImplVertAdv ) THEN
660              CALL GAD_IMPLICIT_R(
661         I         tempImplVertAdv, tempAdvScheme, GAD_TEMPERATURE,
662         I         kappaRk, wVel, theta,
663         U         gT,
664         I         bi, bj, myTime, myIter, myThid )
665            ELSEIF ( tempStepping .AND. implicitDiffusion ) THEN
666    #else /* INCLUDE_IMPLVERTADV_CODE */
667            IF     ( tempStepping .AND. implicitDiffusion ) THEN
668    #endif /* INCLUDE_IMPLVERTADV_CODE */
669  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
670              idkey = iikey + 1  CADJ STORE kappaRk(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
671  CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
672  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
673              CALL IMPLDIFF(            CALL IMPLDIFF(
674       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
675       I         deltaTtracer, KappaRT, recip_HFacC,       I         GAD_TEMPERATURE, kappaRk, recip_hFacC,
676       U         gT,       U         gT,
677       I         myThid )       I         myThid )
678            ENDIF
679    
680    #ifdef ALLOW_TIMEAVE
681            useVariableK = useKPP .OR. usePP81 .OR. useMY82 .OR. useGGL90
682         &       .OR. useGMredi .OR. ivdc_kappa.NE.0.
683            IF (taveFreq.GT.0. .AND. useVariableK ) THEN
684             IF (implicitDiffusion) THEN
685               CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, gT, kappaRk,
686         I                        Nr, 3, deltaTclock, bi, bj, myThid)
687    c        ELSE
688    c          CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, theta, kappaRT,
689    c    I                        Nr, 3, deltaTclock, bi, bj, myThid)
690           ENDIF           ENDIF
691            ENDIF
692    #endif /* ALLOW_TIMEAVE */
693    
694            IF ( saltStepping .AND. implicitDiffusion ) THEN
695              CALL CALC_3D_DIFFUSIVITY(
696         I         bi,bj,iMin,iMax,jMin,jMax,
697         I         GAD_SALINITY, useGMredi, useKPP,
698         O         kappaRk,
699         I         myThid)
700            ENDIF
701    
702           IF (saltStepping) THEN  #ifdef INCLUDE_IMPLVERTADV_CODE
703            IF ( saltImplVertAdv ) THEN
704              CALL GAD_IMPLICIT_R(
705         I         saltImplVertAdv, saltAdvScheme, GAD_SALINITY,
706         I         kappaRk, wVel, salt,
707         U         gS,
708         I         bi, bj, myTime, myIter, myThid )
709            ELSEIF ( saltStepping .AND. implicitDiffusion ) THEN
710    #else /* INCLUDE_IMPLVERTADV_CODE */
711            IF     ( saltStepping .AND. implicitDiffusion ) THEN
712    #endif /* INCLUDE_IMPLVERTADV_CODE */
713  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
714           idkey = iikey + 2  CADJ STORE kappaRk(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
715  CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
716  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
717              CALL IMPLDIFF(            CALL IMPLDIFF(
718       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
719       I         deltaTtracer, KappaRS, recip_HFacC,       I         GAD_SALINITY, kappaRk, recip_hFacC,
720       U         gS,       U         gS,
721       I         myThid )       I         myThid )
722           ENDIF          ENDIF
723    
724  #ifdef ALLOW_PASSIVE_TRACER  #ifdef ALLOW_PTRACERS
725           IF (tr1Stepping) THEN          IF     ( usePTRACERS ) THEN
726  #ifdef ALLOW_AUTODIFF_TAMC  C--     Vertical advection/diffusion (implicit) for passive tracers
727  CADJ STORE gTr1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte             CALL PTRACERS_IMPLICIT(
728  #endif /* ALLOW_AUTODIFF_TAMC */       U                             kappaRk,
729            CALL IMPLDIFF(       I                             bi, bj, myTime, myIter, myThid )
730       I      bi, bj, iMin, iMax, jMin, jMax,          ENDIF
731       I      deltaTtracer, KappaRT, recip_HFacC,  #endif /* ALLOW_PTRACERS */
      U      gTr1,  
      I      myThid )  
          ENDIF  
 #endif  
732    
733  #ifdef   ALLOW_OBCS  #ifdef   ALLOW_OBCS
734  C--      Apply open boundary conditions  C--      Apply open boundary conditions
735           IF (useOBCS) THEN          IF ( ( implicitDiffusion
736         &    .OR. tempImplVertAdv
737         &    .OR. saltImplVertAdv
738         &       ) .AND. useOBCS     ) THEN
739             DO K=1,Nr             DO K=1,Nr
740               CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )               CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )
741             ENDDO             ENDDO
742           END IF          ENDIF
743  #endif   /* ALLOW_OBCS */  #endif   /* ALLOW_OBCS */
744    
745  C--     End If implicitDiffusion  #endif /* SINGLE_LAYER_MODE */
         ENDIF  
746    
747  Ccs-  C--   end bi,bj loops.
748         ENDDO         ENDDO
749        ENDDO        ENDDO
750    
751  #ifdef ALLOW_AIM  #ifdef ALLOW_DEBUG
752        IF ( useAIM ) THEN        If (debugMode) THEN
753         CALL AIM_AIM2DYN_EXCHANGES( myTime, myThid )         CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (THERMODYNAMICS)',myThid)
754        ENDIF         CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (THERMODYNAMICS)',myThid)
755         _EXCH_XYZ_R8(gT,myThid)         CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (THERMODYNAMICS)',myThid)
756         _EXCH_XYZ_R8(gS,myThid)         CALL DEBUG_STATS_RL(Nr,theta,'Theta (THERMODYNAMICS)',myThid)
757  #else         CALL DEBUG_STATS_RL(Nr,salt,'Salt (THERMODYNAMICS)',myThid)
758        IF (staggerTimeStep.AND.useCubedSphereExchange) THEN         CALL DEBUG_STATS_RL(Nr,gT,'Gt (THERMODYNAMICS)',myThid)
759         _EXCH_XYZ_R8(gT,myThid)         CALL DEBUG_STATS_RL(Nr,gS,'Gs (THERMODYNAMICS)',myThid)
760         _EXCH_XYZ_R8(gS,myThid)  #ifndef ALLOW_ADAMSBASHFORTH_3
761           CALL DEBUG_STATS_RL(Nr,gtNm1,'GtNm1 (THERMODYNAMICS)',myThid)
762           CALL DEBUG_STATS_RL(Nr,gsNm1,'GsNm1 (THERMODYNAMICS)',myThid)
763    #endif
764    #ifdef ALLOW_PTRACERS
765           IF ( usePTRACERS ) THEN
766             CALL PTRACERS_DEBUG(myThid)
767           ENDIF
768    #endif /* ALLOW_PTRACERS */
769        ENDIF        ENDIF
770  #endif /* ALLOW_AIM */  #endif /* ALLOW_DEBUG */
771    
772    #ifdef ALLOW_DEBUG
773             IF ( debugLevel .GE. debLevB )
774         &    CALL DEBUG_LEAVE('THERMODYNAMICS',myThid)
775    #endif
776    
777    #endif /* ALLOW_GENERIC_ADVDIFF */
778    
779        RETURN        RETURN
780        END        END
Legend:
Removed from v.1.14  
changed lines
  Added in v.1.112
  ViewVC Help
Powered by ViewVC 1.1.22