/[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.18 by adcroft, Tue Mar 5 14:15:34 2002 UTC revision 1.123 by jmc, Sun Aug 17 02:08:24 2008 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    
4    #include "PACKAGES_CONFIG.h"
5  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
6    #ifdef ALLOW_GENERIC_ADVDIFF
7    # include "GAD_OPTIONS.h"
8    #endif
9    
10    #ifdef ALLOW_AUTODIFF_TAMC
11    # ifdef ALLOW_GMREDI
12    #  include "GMREDI_OPTIONS.h"
13    # endif
14    # ifdef ALLOW_KPP
15    #  include "KPP_OPTIONS.h"
16    # endif
17    #endif /* ALLOW_AUTODIFF_TAMC */
18    
19  CBOP  CBOP
20  C     !ROUTINE: THERMODYNAMICS  C     !ROUTINE: THERMODYNAMICS
# Line 9  C     !INTERFACE: Line 22  C     !INTERFACE:
22        SUBROUTINE THERMODYNAMICS(myTime, myIter, myThid)        SUBROUTINE THERMODYNAMICS(myTime, myIter, myThid)
23  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
24  C     *==========================================================*  C     *==========================================================*
25  C     | SUBROUTINE THERMODYNAMICS                                  C     | SUBROUTINE THERMODYNAMICS
26  C     | o Controlling routine for the prognostic part of the        C     | o Controlling routine for the prognostic part of the
27  C     |   thermo-dynamics.                                          C     |   thermo-dynamics.
28  C     *===========================================================  C     *===========================================================
29  C     | The algorithm...  C     | The algorithm...
30  C     |  C     |
# Line 65  C     == Global variables === Line 78  C     == Global variables ===
78  #include "SIZE.h"  #include "SIZE.h"
79  #include "EEPARAMS.h"  #include "EEPARAMS.h"
80  #include "PARAMS.h"  #include "PARAMS.h"
81    #include "RESTART.h"
82  #include "DYNVARS.h"  #include "DYNVARS.h"
83  #include "GRID.h"  #include "GRID.h"
84  #include "GAD.h"  #ifdef ALLOW_GENERIC_ADVDIFF
85  #ifdef ALLOW_PASSIVE_TRACER  # include "GAD.h"
86  #include "TR1.h"  # include "GAD_SOM_VARS.h"
87    #endif
88    #ifdef ALLOW_PTRACERS
89    # include "PTRACERS_SIZE.h"
90    # include "PTRACERS_PARAMS.h"
91    # include "PTRACERS_FIELDS.h"
92  #endif  #endif
93    #ifdef ALLOW_TIMEAVE
94    # include "TIMEAVE_STATV.h"
95    #endif
96    
97  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
98  # include "tamc.h"  # include "tamc.h"
99  # include "tamc_keys.h"  # include "tamc_keys.h"
100  # include "FFIELDS.h"  # include "FFIELDS.h"
101    # include "SURFACE.h"
102    # include "EOS.h"
103  # ifdef ALLOW_KPP  # ifdef ALLOW_KPP
104  #  include "KPP.h"  #  include "KPP.h"
105  # endif  # endif
106  # ifdef ALLOW_GMREDI  # ifdef ALLOW_GMREDI
107  #  include "GMREDI.h"  #  include "GMREDI.h"
108  # endif  # endif
109    # ifdef ALLOW_EBM
110    #  include "EBM.h"
111    # endif
112  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
 #ifdef ALLOW_TIMEAVE  
 #include "TIMEAVE_STATV.h"  
 #endif  
113    
114  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
115  C     == Routine arguments ==  C     == Routine arguments ==
# Line 95  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
127  C     uTrans, vTrans, rTrans - Per block temporaries holding flow  C     uFld, vFld, wFld       - Local copy of velocity field (3 components)
128  C                              transport  C     uTrans, vTrans, rTrans - Per block temporaries holding flow 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
147  C     k, kup,        - Index for layer above and below. kup and kDown  C     k, kup,        - Index for layer above and below. kup and kDown
148  C     kDown, km1       are switched with layer to be the appropriate  C     kDown, km1       are switched with layer to be the appropriate
149  C                      index into fVerTerm.  C                      index into fVerTerm.
150        _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
151        _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
152          _RL uFld    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
153          _RL vFld    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
154          _RL wFld    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
155        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
156        _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
157        _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
158          _RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
159        _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
160        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
161        _RL fVerS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
162        _RL fVerTr1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL kappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
163        _RL phiHyd  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL kappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
164        _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  #ifdef ALLOW_PTRACERS
165        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fVerP   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2,PTRACERS_num)
166        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL kappaRTr(1-Olx:sNx+Olx,1-Oly:sNy+Oly,PTRACERS_num)
167        _RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  #endif
168        _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)  
169        INTEGER iMin, iMax        INTEGER iMin, iMax
170        INTEGER jMin, jMax        INTEGER jMin, jMax
171        INTEGER bi, bj        INTEGER bi, bj
172        INTEGER i, j        INTEGER i, j
173        INTEGER k, km1, kup, kDown        INTEGER k, km1, kup, kDown
174    #ifdef ALLOW_ADAMSBASHFORTH_3
175          INTEGER iterNb, m1, m2
176    #endif
177    #ifdef ALLOW_TIMEAVE
178          LOGICAL useVariableK
179    #endif
180    #ifdef ALLOW_PTRACERS
181          INTEGER iTracer, ip
182    #endif
183    
184  CEOP  CEOP
185    
186    #ifdef ALLOW_DEBUG
187             IF ( debugLevel .GE. debLevB )
188         &    CALL DEBUG_ENTER('THERMODYNAMICS',myThid)
189    #endif
190    
191  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
192  C--   dummy statement to end declaration part  C--   dummy statement to end declaration part
193        ikey = 1        ikey = 1
194          itdkey = 1
195  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
196    
 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  
   
   
197  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
198  C--   HPF directive to help TAMC  C--   HPF directive to help TAMC
199  CHPF$ INDEPENDENT  CHPF$ INDEPENDENT
200  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
201    
202    C-- Compute correction at the surface for Lin Free Surf.
203    #ifdef ALLOW_AUTODIFF_TAMC
204          TsurfCor = 0. _d 0
205          SsurfCor = 0. _d 0
206    #endif
207          IF (linFSConserveTr) THEN
208    #ifdef ALLOW_AUTODIFF_TAMC
209    CADJ STORE theta,salt,wvel = comlev1, key = ikey_dynamics
210    #endif
211           CALL CALC_WSURF_TR(theta,salt,wVel,
212         &                    myTime,myIter,myThid)
213          ENDIF
214    
215        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
216    
217  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
218  C--    HPF directive to help TAMC  C--    HPF directive to help TAMC
219  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS
220  CHPF$&                  ,phiHyd,utrans,vtrans,xA,yA  CHPF$&                  ,utrans,vtrans,xA,yA
221  CHPF$&                  ,KappaRT,KappaRS  CHPF$&                  ,kappaRT,kappaRS
222  CHPF$&                  )  CHPF$&                  )
223    # ifdef ALLOW_PTRACERS
224    CHPF$  INDEPENDENT, NEW (fVerP,kappaRTr)
225    # endif
226  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
227    
228         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
# Line 208  CHPF$&                  ) Line 235  CHPF$&                  )
235            act3 = myThid - 1            act3 = myThid - 1
236            max3 = nTx*nTy            max3 = nTx*nTy
237            act4 = ikey_dynamics - 1            act4 = ikey_dynamics - 1
238            ikey = (act1 + 1) + act2*max1            itdkey = (act1 + 1) + act2*max1
239       &                      + act3*max1*max2       &                      + act3*max1*max2
240       &                      + act4*max1*max2*max3       &                      + act4*max1*max2*max3
241  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
242    
243  C--     Set up work arrays that need valid initial values  C--   Set up work arrays with valid (i.e. not NaN) values
244    C     These inital values do not alter the numerical results. They
245    C     just ensure that all memory references are to valid floating
246    C     point numbers. This prevents spurious hardware signals due to
247    C     uninitialised but inert locations.
248    
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              xA(i,j)        = 0. _d 0
252              yA(i,j)        = 0. _d 0
253              uTrans(i,j)    = 0. _d 0
254              vTrans(i,j)    = 0. _d 0
255            rTrans (i,j)   = 0. _d 0            rTrans (i,j)   = 0. _d 0
256              rTransKp1(i,j) = 0. _d 0
257            fVerT  (i,j,1) = 0. _d 0            fVerT  (i,j,1) = 0. _d 0
258            fVerT  (i,j,2) = 0. _d 0            fVerT  (i,j,2) = 0. _d 0
259            fVerS  (i,j,1) = 0. _d 0            fVerS  (i,j,1) = 0. _d 0
260            fVerS  (i,j,2) = 0. _d 0            fVerS  (i,j,2) = 0. _d 0
261            fVerTr1(i,j,1) = 0. _d 0            kappaRT(i,j)   = 0. _d 0
262            fVerTr1(i,j,2) = 0. _d 0            kappaRS(i,j)   = 0. _d 0
263           ENDDO           ENDDO
264          ENDDO          ENDDO
265    
# Line 230  C--     Set up work arrays that need val Line 267  C--     Set up work arrays that need val
267           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
268            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
269  C This is currently also used by IVDC and Diagnostics  C This is currently also used by IVDC and Diagnostics
270             ConvectCount(i,j,k) = 0.             kappaRk(i,j,k)    = 0. _d 0
271             KappaRT(i,j,k) = 0. _d 0  C- tracer tendency needs to be set to zero (moved here from gad_calc_rhs):
272             KappaRS(i,j,k) = 0. _d 0             gT(i,j,k,bi,bj)   = 0. _d 0
273  #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  
274            ENDDO            ENDDO
275           ENDDO           ENDDO
276          ENDDO          ENDDO
277    
278          iMin = 1-OLx+1  #ifdef ALLOW_PTRACERS
279          iMax = sNx+OLx          IF ( usePTRACERS ) THEN
280          jMin = 1-OLy+1           DO ip=1,PTRACERS_num
281          jMax = sNy+OLy             DO j=1-OLy,sNy+OLy
282                DO i=1-OLx,sNx+OLx
283                 fVerP  (i,j,1,ip) = 0. _d 0
284  #ifdef ALLOW_AUTODIFF_TAMC               fVerP  (i,j,2,ip) = 0. _d 0
285  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte               kappaRTr(i,j,ip)  = 0. _d 0
286  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte              ENDDO
287  #ifdef ALLOW_KPP             ENDDO
288  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte           ENDDO
289  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  C-      set tracer tendency to zero:
290             DO iTracer=1,PTRACERS_num
291              DO k=1,Nr
292               DO j=1-OLy,sNy+OLy
293                DO i=1-OLx,sNx+OLx
294                 gPTr(i,j,k,bi,bj,itracer) = 0. _d 0
295                ENDDO
296               ENDDO
297              ENDDO
298             ENDDO
299            ENDIF
300  #endif  #endif
 #endif /* ALLOW_AUTODIFF_TAMC */  
301    
302  C--     Start of diagnostic loop  #ifdef ALLOW_ADAMSBASHFORTH_3
303          DO k=Nr,1,-1  C-      Apply AB on T,S :
304            iterNb = myIter
305  #ifdef ALLOW_AUTODIFF_TAMC          IF (staggerTimeStep) iterNb = myIter - 1
306  C? Patrick, is this formula correct now that we change the loop range?          m1 = 1 + MOD(iterNb+1,2)
307  C? Do we still need this?          m2 = 1 + MOD( iterNb ,2)
308  cph kkey formula corrected.  C       compute T^n+1/2 (stored in gtNm) extrapolating T forward in time
309  cph Needed for rhok, rhokm1, in the case useGMREDI.          IF ( AdamsBashforth_T ) CALL ADAMS_BASHFORTH3(
310           kkey = (ikey-1)*Nr + k       I                                  bi, bj, 0,
311  CADJ STORE rhokm1(:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte       U                                  theta, gtNm,
312  CADJ STORE rhok  (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte       I                                  tempStartAB, iterNb, myThid )
313  #endif /* ALLOW_AUTODIFF_TAMC */  C       compute S^n+1/2 (stored in gsNm) extrapolating S forward in time
314            IF ( AdamsBashforth_S ) CALL ADAMS_BASHFORTH3(
315  C--       Integrate continuity vertically for vertical velocity       I                                  bi, bj, 0,
316            CALL INTEGRATE_FOR_W(       U                                  salt, gsNm,
317       I                         bi, bj, k, uVel, vVel,       I                                  saltStartAB, iterNb, myThid )
318       O                         wVel,  #endif /* ALLOW_ADAMSBASHFORTH_3 */
319       I                         myThid )  
320    c       iMin = 1-OLx
321  #ifdef    ALLOW_OBCS  c       iMax = sNx+OLx
322  #ifdef    ALLOW_NONHYDROSTATIC  c       jMin = 1-OLy
323  C--       Apply OBC to W if in N-H mode  c       jMax = sNy+OLy
           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  
324    
325  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
326  cph avoids recomputation of integrate_for_w  cph avoids recomputation of integrate_for_w
327  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  
328  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
329    
330  #ifdef  ALLOW_GMREDI  C--     Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h
331    C--     MOST of THERMODYNAMICS will be disabled
332    #ifndef SINGLE_LAYER_MODE
333    
334  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
335  CADJ STORE sigmaX(:,:,:) = comlev1, key=ikey, byte=isbyte  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
336  CADJ STORE sigmaY(:,:,:) = comlev1, key=ikey, byte=isbyte  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
337  CADJ STORE sigmaR(:,:,:) = comlev1, key=ikey, byte=isbyte  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
338  #endif /* ALLOW_AUTODIFF_TAMC */  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
339  C--     Calculate iso-neutral slopes for the GM/Redi parameterisation  # if ((defined ALLOW_DEPTH_CONTROL) || (defined NONLIN_FRSURF))
340          IF (useGMRedi) THEN  CADJ STORE gtnm1(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
341            CALL GMREDI_CALC_TENSOR(  CADJ STORE gsnm1(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
342       I             bi, bj, iMin, iMax, jMin, jMax,  # endif
      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 */  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE KappaRT(:,:,:)     = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE KappaRS(:,:,:)     = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  
 #ifdef ALLOW_PASSIVE_TRACER  
 CADJ STORE tr1  (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  
 #endif  
343  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
344    
 #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 */  
   
345  #ifndef DISABLE_MULTIDIM_ADVECTION  #ifndef DISABLE_MULTIDIM_ADVECTION
346  C--     Some advection schemes are better calculated using a multi-dimensional  C--     Some advection schemes are better calculated using a multi-dimensional
347  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 454  C to be able to exclude this scheme to a Line 353  C to be able to exclude this scheme to a
353  C recomputation. It *is* differentiable, if you need it.  C recomputation. It *is* differentiable, if you need it.
354  C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to  C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to
355  C disable this section of code.  C disable this section of code.
356          IF (multiDimAdvection) THEN  #ifdef GAD_ALLOW_TS_SOM_ADV
357           IF (tempStepping .AND.          IF ( tempSOM_Advection ) THEN
358       &       tempAdvScheme.NE.ENUM_CENTERED_2ND .AND.  #ifdef ALLOW_DEBUG
359       &       tempAdvScheme.NE.ENUM_UPWIND_3RD .AND.            IF ( debugLevel .GE. debLevB )
360       &       tempAdvScheme.NE.ENUM_CENTERED_4TH ) THEN       &     CALL DEBUG_CALL('GAD_SOM_ADVECT',myThid)
361            CALL GAD_ADVECTION(bi,bj,tempAdvScheme,GAD_TEMPERATURE,  #endif
362       U                      theta,gT,            CALL GAD_SOM_ADVECT(
363       I                      myTime,myIter,myThid)       I             tempImplVertAdv, tempAdvScheme, tempVertAdvScheme,
364           ENDIF       I             GAD_TEMPERATURE,
365           IF (saltStepping .AND.       I             uVel, vVel, wVel, theta,
366       &       saltAdvScheme.NE.ENUM_CENTERED_2ND .AND.       U             som_T,
367       &       saltAdvScheme.NE.ENUM_UPWIND_3RD .AND.       O             gT,
368       &       saltAdvScheme.NE.ENUM_CENTERED_4TH ) THEN       I             bi,bj,myTime,myIter,myThid)
369            CALL GAD_ADVECTION(bi,bj,saltAdvScheme,GAD_SALINITY,          ELSEIF (tempMultiDimAdvec) THEN
370       U                      salt,gS,  #else /* GAD_ALLOW_TS_SOM_ADV */
371       I                      myTime,myIter,myThid)          IF (tempMultiDimAdvec) THEN
372           ENDIF  #endif /* GAD_ALLOW_TS_SOM_ADV */
373    #ifdef ALLOW_DEBUG
374              IF ( debugLevel .GE. debLevB )
375         &     CALL DEBUG_CALL('GAD_ADVECTION',myThid)
376    #endif
377              CALL GAD_ADVECTION(
378         I             tempImplVertAdv, tempAdvScheme, tempVertAdvScheme,
379         I             GAD_TEMPERATURE,
380         I             uVel, vVel, wVel, theta,
381         O             gT,
382         I             bi,bj,myTime,myIter,myThid)
383          ENDIF          ENDIF
384    #ifdef GAD_ALLOW_TS_SOM_ADV
385            IF ( saltSOM_Advection ) THEN
386    #ifdef ALLOW_DEBUG
387              IF ( debugLevel .GE. debLevB )
388         &     CALL DEBUG_CALL('GAD_SOM_ADVECT',myThid)
389    #endif
390              CALL GAD_SOM_ADVECT(
391         I             saltImplVertAdv, saltAdvScheme, saltVertAdvScheme,
392         I             GAD_SALINITY,
393         I             uVel, vVel, wVel, salt,
394         U             som_S,
395         O             gS,
396         I             bi,bj,myTime,myIter,myThid)
397            ELSEIF (saltMultiDimAdvec) THEN
398    #else /* GAD_ALLOW_TS_SOM_ADV */
399            IF (saltMultiDimAdvec) THEN
400    #endif /* GAD_ALLOW_TS_SOM_ADV */
401    #ifdef ALLOW_DEBUG
402              IF ( debugLevel .GE. debLevB )
403         &     CALL DEBUG_CALL('GAD_ADVECTION',myThid)
404    #endif
405              CALL GAD_ADVECTION(
406         I             saltImplVertAdv, saltAdvScheme, saltVertAdvScheme,
407         I             GAD_SALINITY,
408         I             uVel, vVel, wVel, salt,
409         O             gS,
410         I             bi,bj,myTime,myIter,myThid)
411            ENDIF
412    
413  C Since passive tracers are configurable separately from T,S we  C Since passive tracers are configurable separately from T,S we
414  C call the multi-dimensional method for PTRACERS regardless  C call the multi-dimensional method for PTRACERS regardless
415  C of whether multiDimAdvection is set or not.  C of whether multiDimAdvection is set or not.
416  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
417          IF ( usePTRACERS ) THEN          IF ( usePTRACERS ) THEN
418    #ifdef ALLOW_DEBUG
419              IF ( debugLevel .GE. debLevB )
420         &     CALL DEBUG_CALL('PTRACERS_ADVECTION',myThid)
421    #endif
422           CALL PTRACERS_ADVECTION( bi,bj,myIter,myTime,myThid )           CALL PTRACERS_ADVECTION( bi,bj,myIter,myTime,myThid )
423          ENDIF          ENDIF
424  #endif /* ALLOW_PTRACERS */  #endif /* ALLOW_PTRACERS */
425  #endif /* DISABLE_MULTIDIM_ADVECTION */  #endif /* DISABLE_MULTIDIM_ADVECTION */
426    
427    #ifdef ALLOW_DEBUG
428            IF ( debugLevel .GE. debLevB )
429         &    CALL DEBUG_MSG('ENTERING DOWNWARD K LOOP',myThid)
430    #endif
431    
432  C--     Start of thermodynamics loop  C--     Start of thermodynamics loop
433          DO k=Nr,1,-1          DO k=Nr,1,-1
434  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
435  C? Patrick Is this formula correct?  C? Patrick Is this formula correct?
436  cph Yes, but I rewrote it.  cph Yes, but I rewrote it.
437  cph Also, the KappaR? need the index and subscript k!  cph Also, the kappaR? need the index and subscript k!
438           kkey = (ikey-1)*Nr + k           kkey = (itdkey-1)*Nr + k
439  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
440    
441  C--       km1    Points to level above k (=k-1)  C--       km1    Points to level above k (=k-1)
# Line 504  C--       kDown  Cycles through 2,1 to p Line 451  C--       kDown  Cycles through 2,1 to p
451            jMin = 1-OLy            jMin = 1-OLy
452            jMax = sNy+OLy            jMax = sNy+OLy
453    
454  C--       Get temporary terms used by tendency routines            IF (k.EQ.Nr) THEN
455               DO j=1-Oly,sNy+Oly
456                DO i=1-Olx,sNx+Olx
457                 rTransKp1(i,j) = 0. _d 0
458                ENDDO
459               ENDDO
460              ELSE
461               DO j=1-Oly,sNy+Oly
462                DO i=1-Olx,sNx+Olx
463                 rTransKp1(i,j) = rTrans(i,j)
464                ENDDO
465               ENDDO
466              ENDIF
467    #ifdef ALLOW_AUTODIFF_TAMC
468    CADJ STORE rTransKp1(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
469    #endif
470    
471    C--       Get temporary terms used by tendency routines :
472    C-        Calculate horizontal "volume transport" through tracer cell face
473    C         anelastic: uTrans,vTrans are scaled by rhoFacC (~ mass transport)
474            CALL CALC_COMMON_FACTORS (            CALL CALC_COMMON_FACTORS (
475       I         bi,bj,iMin,iMax,jMin,jMax,k,       I         uVel, vVel,
476       O         xA,yA,uTrans,vTrans,rTrans,maskUp,       O         uFld, vFld, uTrans, vTrans, xA, yA,
477       I         myThid)       I         k,bi,bj, myThid )
478    
479    C-        Calculate vertical "volume transport" through tracer cell face
480              IF (k.EQ.1) THEN
481    C-         Surface interface :
482               DO j=1-Oly,sNy+Oly
483                DO i=1-Olx,sNx+Olx
484                 wFld(i,j)   = 0. _d 0
485                 maskUp(i,j) = 0. _d 0
486                 rTrans(i,j) = 0. _d 0
487                ENDDO
488               ENDDO
489              ELSE
490    C-         Interior interface :
491    C          anelastic: rTrans is scaled by rhoFacF (~ mass transport)
492               DO j=1-Oly,sNy+Oly
493                DO i=1-Olx,sNx+Olx
494                 wFld(i,j)   = wVel(i,j,k,bi,bj)
495                 maskUp(i,j) = maskC(i,j,k-1,bi,bj)*maskC(i,j,k,bi,bj)
496                 rTrans(i,j) = wFld(i,j)*rA(i,j,bi,bj)*maskUp(i,j)
497         &                              *deepFac2F(k)*rhoFacF(k)
498                ENDDO
499               ENDDO
500              ENDIF
501    
502  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_GMREDI
503  CADJ STORE KappaRT(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte  C--   Residual transp = Bolus transp + Eulerian transp
504  CADJ STORE KappaRS(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte            IF (useGMRedi) THEN
505  #endif /* ALLOW_AUTODIFF_TAMC */              CALL GMREDI_CALC_UVFLOW(
506         U                  uFld, vFld, uTrans, vTrans,
507         I                  k, bi, bj, myThid )
508                IF (K.GE.2) THEN
509                  CALL GMREDI_CALC_WFLOW(
510         U                  wFld, rTrans,
511         I                  k, bi, bj, myThid )
512                ENDIF
513              ENDIF
514    # ifdef ALLOW_AUTODIFF_TAMC
515    CADJ STORE rTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
516    CADJ STORE wfld(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
517    # ifdef GM_BOLUS_ADVEC
518    CADJ STORE ufld(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
519    CADJ STORE vfld(:,:)      = comlev1_bibj_k, key=kkey, byte=isbyte
520    CADJ STORE uTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
521    CADJ STORE vTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
522    # endif
523    # endif /* ALLOW_AUTODIFF_TAMC */
524    #endif /* ALLOW_GMREDI */
525    
526  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL
527  C--      Calculate the total vertical diffusivity  C--      Calculate the total vertical diffusivity
528           CALL CALC_DIFFUSIVITY(            IF ( .NOT.implicitDiffusion ) THEN
529       I        bi,bj,iMin,iMax,jMin,jMax,k,              CALL CALC_DIFFUSIVITY(
530       I        maskUp,       I          bi,bj,iMin,iMax,jMin,jMax,k,
531       O        KappaRT,KappaRS,       I          maskUp,
532       I        myThid)       O          kappaRT,kappaRS,
533         I          myThid)
534              ENDIF
535    # ifdef ALLOW_AUTODIFF_TAMC
536    CADJ STORE kappaRT(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
537    CADJ STORE kappaRS(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
538    # endif /* ALLOW_AUTODIFF_TAMC */
539  #endif  #endif
540    
541            iMin = 1-OLx+2            iMin = 1-OLx+2
# Line 530  C--      Calculate the total vertical di Line 544  C--      Calculate the total vertical di
544            jMax = sNy+OLy-1            jMax = sNy+OLy-1
545    
546  C--      Calculate active tracer tendencies (gT,gS,...)  C--      Calculate active tracer tendencies (gT,gS,...)
547  C        and step forward storing result in gTnm1, gSnm1, etc.  C        and step forward storing result in gT, gS, etc.
548    C--
549    # ifdef ALLOW_AUTODIFF_TAMC
550    #  if ((defined NONLIN_FRSURF) || (defined ALLOW_DEPTH_CONTROL)) && (defined ALLOW_GMREDI)
551    #   ifdef GM_NON_UNITY_DIAGONAL
552    CADJ STORE kux(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
553    CADJ STORE kvy(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
554    #   endif
555    #   ifdef GM_EXTRA_DIAGONAL
556    CADJ STORE kuz(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
557    CADJ STORE kvz(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
558    #   endif
559    #  endif
560    # endif /* ALLOW_AUTODIFF_TAMC */
561    C
562    #ifdef ALLOW_AUTODIFF_TAMC
563    # if (defined NONLIN_FRSURF) || (defined ALLOW_DEPTH_CONTROL)
564    cph-test
565    CADJ STORE uFld(:,:), vFld(:,:), wFld(:,:)
566    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
567    CADJ STORE uTrans(:,:), vTrans(:,:)
568    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
569    CADJ STORE xA(:,:), yA(:,:)
570    CADJ &     = comlev1_bibj_k, key=kkey, byte=isbyte
571    # endif
572    #endif /* ALLOW_AUTODIFF_TAMC */
573    C
574           IF ( tempStepping ) THEN           IF ( tempStepping ) THEN
575             CALL CALC_GT(  #ifdef ALLOW_AUTODIFF_TAMC
576    CADJ STORE gTnm1(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
577    # if (defined NONLIN_FRSURF) || (defined ALLOW_DEPTH_CONTROL)
578    CADJ STORE gt(:,:,:,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
579    CADJ STORE fvert(:,:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
580    # endif
581    #endif /* ALLOW_AUTODIFF_TAMC */
582              CALL CALC_GT(
583       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
584       I         xA,yA,uTrans,vTrans,rTrans,maskUp,       I         xA, yA, maskUp, uFld, vFld, wFld,
585       I         KappaRT,       I         uTrans, vTrans, rTrans, rTransKp1,
586         I         kappaRT,
587       U         fVerT,       U         fVerT,
588       I         myTime,myIter,myThid)       I         myTime,myIter,myThid)
589    #ifdef ALLOW_ADAMSBASHFORTH_3
590              IF ( AdamsBashforth_T ) THEN
591             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
592       I         bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,       I         bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,
593       I         theta, gT,       I         gtNm(1-Olx,1-Oly,1,1,1,m2),
594         U         gT,
595         I         myIter, myThid)
596              ELSE
597    #endif
598               CALL TIMESTEP_TRACER(
599         I         bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,
600         I         theta,
601         U         gT,
602       I         myIter, myThid)       I         myIter, myThid)
603    #ifdef ALLOW_ADAMSBASHFORTH_3
604              ENDIF
605    #endif
606           ENDIF           ENDIF
607    
608           IF ( saltStepping ) THEN           IF ( saltStepping ) THEN
609             CALL CALC_GS(  #ifdef ALLOW_AUTODIFF_TAMC
610    CADJ STORE gSnm1(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
611    # if (defined NONLIN_FRSURF) || (defined ALLOW_DEPTH_CONTROL)
612    CADJ STORE gs(:,:,:,bi,bj)    = comlev1_bibj_k, key=kkey, byte=isbyte
613    CADJ STORE fvers(:,:,:)       = comlev1_bibj_k, key=kkey, byte=isbyte
614    # endif
615    #endif /* ALLOW_AUTODIFF_TAMC */
616              CALL CALC_GS(
617       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
618       I         xA,yA,uTrans,vTrans,rTrans,maskUp,       I         xA, yA, maskUp, uFld, vFld, wFld,
619       I         KappaRS,       I         uTrans, vTrans, rTrans, rTransKp1,
620         I         kappaRS,
621       U         fVerS,       U         fVerS,
622       I         myTime,myIter,myThid)       I         myTime,myIter,myThid)
623    #ifdef ALLOW_ADAMSBASHFORTH_3
624              IF ( AdamsBashforth_S ) THEN
625             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
626       I         bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,       I         bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,
627       I         salt, gS,       I         gsNm(1-Olx,1-Oly,1,1,1,m2),
628         U         gS,
629       I         myIter, myThid)       I         myIter, myThid)
630           ENDIF            ELSE
631  #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)  
632             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
633       I         bi,bj,iMin,iMax,jMin,jMax,k,tracerAdvScheme,       I         bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,
634       I         Tr1, gTr1,       I         salt,
635       I         myIter,myThid)       U         gS,
636           ENDIF       I         myIter, myThid)
637    #ifdef ALLOW_ADAMSBASHFORTH_3
638              ENDIF
639  #endif  #endif
640             ENDIF
641    
642  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
643           IF ( usePTRACERS ) THEN           IF ( usePTRACERS ) THEN
644             CALL PTRACERS_INTEGERATE(             IF ( .NOT.implicitDiffusion ) THEN
645       I         bi,bj,k,               CALL PTRACERS_CALC_DIFF(
646       I         xA,yA,uTrans,vTrans,rTrans,maskUp,       I            bi,bj,iMin,iMax,jMin,jMax,k,
647       X         KappaRS,       I            maskUp,
648       I         myIter,myTime,myThid)       O            kappaRTr,
649         I            myThid)
650               ENDIF
651    # ifdef ALLOW_AUTODIFF_TAMC
652    CADJ STORE kappaRTr(:,:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
653    # endif /* ALLOW_AUTODIFF_TAMC */
654               CALL PTRACERS_INTEGRATE(
655         I          bi,bj,k,
656         I          xA, yA, maskUp, uFld, vFld, wFld,
657         I          uTrans, vTrans, rTrans, rTransKp1,
658         I          kappaRTr,
659         U          fVerP,
660         I          myTime,myIter,myThid)
661           ENDIF           ENDIF
662  #endif /* ALLOW_PTRACERS */  #endif /* ALLOW_PTRACERS */
663    
# Line 587  C--      Apply open boundary conditions Line 669  C--      Apply open boundary conditions
669  #endif   /* ALLOW_OBCS */  #endif   /* ALLOW_OBCS */
670    
671  C--      Freeze water  C--      Freeze water
672           IF (allowFreezing) THEN  C  this bit of code is left here for backward compatibility.
673    C  freezing at surface level has been moved to FORWARD_STEP
674             IF ( useOldFreezing .AND. .NOT. useSEAICE
675         &       .AND. .NOT.(useThSIce.AND.k.EQ.1) ) THEN
676  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
677  CADJ STORE gT(:,:,k,bi,bj) = comlev1_bibj_k  CADJ STORE gT(:,:,k,bi,bj) = comlev1_bibj_k
678  CADJ &   , key = kkey, byte = isbyte  CADJ &   , key = kkey, byte = isbyte
679  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
680              CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )              CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )
681           END IF           ENDIF
682    
683  C--     end of thermodynamic k loop (Nr:1)  C--     end of thermodynamic k loop (Nr:1)
684          ENDDO          ENDDO
685    
686    #ifdef ALLOW_DOWN_SLOPE
687            IF ( tempStepping .AND. useDOWN_SLOPE ) THEN
688              IF ( usingPCoords ) THEN
689                CALL DWNSLP_APPLY(
690         I                  GAD_TEMPERATURE, bi, bj, kSurfC,
691         I                  recip_drF, recip_hFacC, recip_rA,
692         I                  theta,
693         U                  gT,
694         I                  myTime, myIter, myThid )
695              ELSE
696                CALL DWNSLP_APPLY(
697         I                  GAD_TEMPERATURE, bi, bj, kLowC,
698         I                  recip_drF, recip_hFacC, recip_rA,
699         I                  theta,
700         U                  gT,
701         I                  myTime, myIter, myThid )
702              ENDIF
703            ENDIF
704            IF ( saltStepping .AND. useDOWN_SLOPE ) THEN
705              IF ( usingPCoords ) THEN
706                CALL DWNSLP_APPLY(
707         I                  GAD_SALINITY, bi, bj, kSurfC,
708         I                  recip_drF, recip_hFacC, recip_rA,
709         I                  salt,
710         U                  gS,
711         I                  myTime, myIter, myThid )
712              ELSE
713                CALL DWNSLP_APPLY(
714         I                  GAD_SALINITY, bi, bj, kLowC,
715         I                  recip_drF, recip_hFacC, recip_rA,
716         I                  salt,
717         U                  gS,
718         I                  myTime, myIter, myThid )
719              ENDIF
720            ENDIF
721    #ifdef ALLOW_PTRACERS
722            IF ( usePTRACERS .AND. useDOWN_SLOPE ) THEN
723              CALL PTRACERS_DWNSLP_APPLY(
724         I                  bi, bj, myTime, myIter, myThid )
725            ENDIF
726    #endif /* ALLOW_PTRACERS */
727    #endif /* ALLOW_DOWN_SLOPE */
728    
729  #ifdef ALLOW_AUTODIFF_TAMC  C       All explicit advection/diffusion/sources should now be
730  C? Patrick? What about this one?  C       done. The updated tracer field is in gPtr. Accumalate
731  cph Keys iikey and idkey dont seem to be needed  C       explicit tendency and also reset gPtr to initial tracer
732  cph since storing occurs on different tape for each  C       field for implicit matrix calculation
733  cph impldiff call anyways.  
734  cph Thus, common block comlev1_impl isnt needed either.  #ifdef ALLOW_MATRIX
735  cph Storing below needed in the case useGMREDI.          IF (useMATRIX)
736          iikey = (ikey-1)*maximpl       &    CALL MATRIX_STORE_TENDENCY_EXP(bi,bj, myTime,myIter,myThid)
737  #endif /* ALLOW_AUTODIFF_TAMC */  #endif
   
 C--     Implicit diffusion  
         IF (implicitDiffusion) THEN  
738    
739           IF (tempStepping) THEN          iMin = 1
740            iMax = sNx
741            jMin = 1
742            jMax = sNy
743    
744    C--     Implicit vertical advection & diffusion
745            IF ( tempStepping .AND. implicitDiffusion ) THEN
746              CALL CALC_3D_DIFFUSIVITY(
747         I         bi,bj,iMin,iMax,jMin,jMax,
748         I         GAD_TEMPERATURE, useGMredi, useKPP,
749         O         kappaRk,
750         I         myThid)
751            ENDIF
752    #ifdef INCLUDE_IMPLVERTADV_CODE
753            IF ( tempImplVertAdv ) THEN
754              CALL GAD_IMPLICIT_R(
755         I         tempImplVertAdv, tempAdvScheme, GAD_TEMPERATURE,
756         I         kappaRk, wVel, theta,
757         U         gT,
758         I         bi, bj, myTime, myIter, myThid )
759            ELSEIF ( tempStepping .AND. implicitDiffusion ) THEN
760    #else /* INCLUDE_IMPLVERTADV_CODE */
761            IF     ( tempStepping .AND. implicitDiffusion ) THEN
762    #endif /* INCLUDE_IMPLVERTADV_CODE */
763  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
764              idkey = iikey + 1  CADJ STORE kappaRk(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
765  CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
766  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
767              CALL IMPLDIFF(            CALL IMPLDIFF(
768       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
769       I         deltaTtracer, KappaRT, recip_HFacC,       I         GAD_TEMPERATURE, kappaRk, recip_hFacC,
770       U         gT,       U         gT,
771       I         myThid )       I         myThid )
772            ENDIF
773    
774    #ifdef ALLOW_TIMEAVE
775            useVariableK = useKPP .OR. usePP81 .OR. useMY82 .OR. useGGL90
776         &       .OR. useGMredi .OR. ivdc_kappa.NE.0.
777            IF (taveFreq.GT.0. .AND. useVariableK ) THEN
778             IF (implicitDiffusion) THEN
779               CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, gT, kappaRk,
780         I                        Nr, 3, deltaTclock, bi, bj, myThid)
781    c        ELSE
782    c          CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, theta, kappaRT,
783    c    I                        Nr, 3, deltaTclock, bi, bj, myThid)
784           ENDIF           ENDIF
785            ENDIF
786    #endif /* ALLOW_TIMEAVE */
787    
788           IF (saltStepping) THEN          IF ( saltStepping .AND. implicitDiffusion ) THEN
789              CALL CALC_3D_DIFFUSIVITY(
790         I         bi,bj,iMin,iMax,jMin,jMax,
791         I         GAD_SALINITY, useGMredi, useKPP,
792         O         kappaRk,
793         I         myThid)
794            ENDIF
795    
796    #ifdef INCLUDE_IMPLVERTADV_CODE
797            IF ( saltImplVertAdv ) THEN
798              CALL GAD_IMPLICIT_R(
799         I         saltImplVertAdv, saltAdvScheme, GAD_SALINITY,
800         I         kappaRk, wVel, salt,
801         U         gS,
802         I         bi, bj, myTime, myIter, myThid )
803            ELSEIF ( saltStepping .AND. implicitDiffusion ) THEN
804    #else /* INCLUDE_IMPLVERTADV_CODE */
805            IF     ( saltStepping .AND. implicitDiffusion ) THEN
806    #endif /* INCLUDE_IMPLVERTADV_CODE */
807  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
808           idkey = iikey + 2  CADJ STORE kappaRk(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
809  CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
810  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
811              CALL IMPLDIFF(            CALL IMPLDIFF(
812       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
813       I         deltaTtracer, KappaRS, recip_HFacC,       I         GAD_SALINITY, kappaRk, recip_hFacC,
814       U         gS,       U         gS,
815       I         myThid )       I         myThid )
816           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  
 #endif  
817    
818  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
819  C Vertical diffusion (implicit) for passive tracers          IF     ( usePTRACERS ) THEN
820           IF ( usePTRACERS ) THEN  C--     Vertical advection/diffusion (implicit) for passive tracers
821             CALL PTRACERS_IMPLDIFF( bi,bj,KappaRS,myThid )             CALL PTRACERS_IMPLICIT(
822           ENDIF       U                             kappaRk,
823         I                             bi, bj, myTime, myIter, myThid )
824            ENDIF
825  #endif /* ALLOW_PTRACERS */  #endif /* ALLOW_PTRACERS */
826    
827  #ifdef   ALLOW_OBCS  #ifdef   ALLOW_OBCS
828  C--      Apply open boundary conditions  C--      Apply open boundary conditions
829           IF (useOBCS) THEN          IF ( ( implicitDiffusion
830         &    .OR. tempImplVertAdv
831         &    .OR. saltImplVertAdv
832         &       ) .AND. useOBCS     ) THEN
833             DO K=1,Nr             DO K=1,Nr
834               CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )               CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )
835             ENDDO             ENDDO
836           END IF          ENDIF
837  #endif   /* ALLOW_OBCS */  #endif   /* ALLOW_OBCS */
838    
839  C--     End If implicitDiffusion  #endif /* SINGLE_LAYER_MODE */
         ENDIF  
840    
841  Ccs-  C--   end bi,bj loops.
842         ENDDO         ENDDO
843        ENDDO        ENDDO
844    
845  #ifdef ALLOW_AIM  #ifdef ALLOW_DEBUG
       IF ( useAIM ) THEN  
        CALL AIM_AIM2DYN_EXCHANGES( myTime, myThid )  
       ENDIF  
        _EXCH_XYZ_R8(gT,myThid)  
        _EXCH_XYZ_R8(gS,myThid)  
 #else  
       IF (staggerTimeStep.AND.useCubedSphereExchange) THEN  
        _EXCH_XYZ_R8(gT,myThid)  
        _EXCH_XYZ_R8(gS,myThid)  
       ENDIF  
 #endif /* ALLOW_AIM */  
   
 #ifndef DISABLE_DEBUGMODE  
846        If (debugMode) THEN        If (debugMode) THEN
847         CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (THERMODYNAMICS)',myThid)
848         CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (THERMODYNAMICS)',myThid)
849         CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (THERMODYNAMICS)',myThid)
850         CALL DEBUG_STATS_RL(Nr,theta,'Theta (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,theta,'Theta (THERMODYNAMICS)',myThid)
851         CALL DEBUG_STATS_RL(Nr,salt,'Salt (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,salt,'Salt (THERMODYNAMICS)',myThid)
852         CALL DEBUG_STATS_RL(Nr,Gt,'Gt (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,gT,'Gt (THERMODYNAMICS)',myThid)
853         CALL DEBUG_STATS_RL(Nr,Gs,'Gs (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,gS,'Gs (THERMODYNAMICS)',myThid)
854         CALL DEBUG_STATS_RL(Nr,GtNm1,'GtNm1 (THERMODYNAMICS)',myThid)  #ifndef ALLOW_ADAMSBASHFORTH_3
855         CALL DEBUG_STATS_RL(Nr,GsNm1,'GsNm1 (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,gtNm1,'GtNm1 (THERMODYNAMICS)',myThid)
856           CALL DEBUG_STATS_RL(Nr,gsNm1,'GsNm1 (THERMODYNAMICS)',myThid)
857    #endif
858  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
859         IF ( usePTRACERS ) THEN         IF ( usePTRACERS ) THEN
860           CALL PTRACERS_DEBUG(myThid)           CALL PTRACERS_DEBUG(myThid)
861         ENDIF         ENDIF
862  #endif /* ALLOW_PTRACERS */  #endif /* ALLOW_PTRACERS */
863        ENDIF        ENDIF
864    #endif /* ALLOW_DEBUG */
865    
866    #ifdef ALLOW_DEBUG
867             IF ( debugLevel .GE. debLevB )
868         &    CALL DEBUG_LEAVE('THERMODYNAMICS',myThid)
869  #endif  #endif
870    
871    #endif /* ALLOW_GENERIC_ADVDIFF */
872    
873        RETURN        RETURN
874        END        END
Legend:
Removed from v.1.18  
changed lines
  Added in v.1.123
  ViewVC Help
Powered by ViewVC 1.1.22