/[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.31 by cheisey, Fri Nov 15 19:58:21 2002 UTC revision 1.105 by jmc, Sun Jun 18 23:22:43 2006 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    
4    #include "PACKAGES_CONFIG.h"
5  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
6    
7  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
8  # ifdef ALLOW_GMREDI  # ifdef ALLOW_GMREDI
9  #  include "GMREDI_OPTIONS.h"  #  include "GMREDI_OPTIONS.h"
# Line 9  C $Name$ Line 11  C $Name$
11  # ifdef ALLOW_KPP  # ifdef ALLOW_KPP
12  #  include "KPP_OPTIONS.h"  #  include "KPP_OPTIONS.h"
13  # endif  # endif
 cswdice --- add ----  
 #ifdef ALLOW_TSEAICE  
 #include "ICE.h"  
 #endif  
 cswdice ------  
14  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
15    
16  CBOP  CBOP
# Line 22  C     !INTERFACE: Line 19  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 80  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"
# Line 95  C     == Global variables === Line 99  C     == Global variables ===
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 109  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 140  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        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  
         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
# Line 201  CHPF$ INDEPENDENT Line 197  CHPF$ INDEPENDENT
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, NEW (rTrans,fVerT,fVerS  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS
200  CHPF$&                  ,phiHyd,utrans,vtrans,xA,yA  CHPF$&                  ,utrans,vtrans,xA,yA
201  CHPF$&                  ,KappaRT,KappaRS  CHPF$&                  ,kappaRT,kappaRS
202  CHPF$&                  )  CHPF$&                  )
203    # ifdef ALLOW_PTRACERS
204    CHPF$  INDEPENDENT, NEW (fVerP,kappaRTr)
205    # endif
206  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
207    
208         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
# Line 221  CHPF$&                  ) Line 220  CHPF$&                  )
220       &                      + act4*max1*max2*max3       &                      + act4*max1*max2*max3
221  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
222    
223  C--     Set up work arrays that need valid initial values  C--   Set up work arrays with valid (i.e. not NaN) values
224    C     These inital values do not alter the numerical results. They
225    C     just ensure that all memory references are to valid floating
226    C     point numbers. This prevents spurious hardware signals due to
227    C     uninitialised but inert locations.
228    
229          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
230           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
231              xA(i,j)        = 0. _d 0
232              yA(i,j)        = 0. _d 0
233              uTrans(i,j)    = 0. _d 0
234              vTrans(i,j)    = 0. _d 0
235            rTrans (i,j)   = 0. _d 0            rTrans (i,j)   = 0. _d 0
236              rTransKp1(i,j) = 0. _d 0
237            fVerT  (i,j,1) = 0. _d 0            fVerT  (i,j,1) = 0. _d 0
238            fVerT  (i,j,2) = 0. _d 0            fVerT  (i,j,2) = 0. _d 0
239            fVerS  (i,j,1) = 0. _d 0            fVerS  (i,j,1) = 0. _d 0
240            fVerS  (i,j,2) = 0. _d 0            fVerS  (i,j,2) = 0. _d 0
241            fVerTr1(i,j,1) = 0. _d 0            kappaRT(i,j)   = 0. _d 0
242            fVerTr1(i,j,2) = 0. _d 0            kappaRS(i,j)   = 0. _d 0
           rhoKM1 (i,j)   = 0. _d 0  
243           ENDDO           ENDDO
244          ENDDO          ENDDO
245    
# Line 239  C--     Set up work arrays that need val Line 247  C--     Set up work arrays that need val
247           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
248            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
249  C This is currently also used by IVDC and Diagnostics  C This is currently also used by IVDC and Diagnostics
250             phiHyd(i,j,k) = 0. _d 0             kappaRk(i,j,k)    = 0. _d 0
251             sigmaX(i,j,k) = 0. _d 0  C- tracer tendency needs to be set to zero (moved here from gad_calc_rhs):
            sigmaY(i,j,k) = 0. _d 0  
            sigmaR(i,j,k) = 0. _d 0  
            ConvectCount(i,j,k) = 0.  
            KappaRT(i,j,k)    = 0. _d 0  
            KappaRS(i,j,k)    = 0. _d 0  
 #ifdef ALLOW_AUTODIFF_TAMC  
 cph all the following init. are necessary for TAF  
 cph although some of these are re-initialised later.  
252             gT(i,j,k,bi,bj)   = 0. _d 0             gT(i,j,k,bi,bj)   = 0. _d 0
253             gS(i,j,k,bi,bj)   = 0. _d 0             gS(i,j,k,bi,bj)   = 0. _d 0
 # ifdef ALLOW_PASSIVE_TRACER  
            gTr1(i,j,k,bi,bj) = 0. _d 0  
 # endif  
 # ifdef ALLOW_GMREDI  
            Kwx(i,j,k,bi,bj)  = 0. _d 0  
            Kwy(i,j,k,bi,bj)  = 0. _d 0  
            Kwz(i,j,k,bi,bj)  = 0. _d 0  
 #  ifdef GM_NON_UNITY_DIAGONAL  
            Kux(i,j,k,bi,bj)  = 0. _d 0  
            Kvy(i,j,k,bi,bj)  = 0. _d 0  
 #  endif  
 #  ifdef GM_EXTRA_DIAGONAL  
            Kuz(i,j,k,bi,bj)  = 0. _d 0  
            Kvz(i,j,k,bi,bj)  = 0. _d 0  
 #  endif  
 #  ifdef GM_BOLUS_ADVEC  
            GM_PsiX(i,j,k,bi,bj)  = 0. _d 0  
            GM_PsiY(i,j,k,bi,bj)  = 0. _d 0  
 #  endif  
 # endif /* ALLOW_GMREDI */  
 #endif /* ALLOW_AUTODIFF_TAMC */  
254            ENDDO            ENDDO
255           ENDDO           ENDDO
256          ENDDO          ENDDO
257    
258          iMin = 1-OLx  #ifdef ALLOW_PTRACERS
259          iMax = sNx+OLx          IF ( usePTRACERS ) THEN
260          jMin = 1-OLy           DO ip=1,PTRACERS_num
261          jMax = sNy+OLy             DO j=1-OLy,sNy+OLy
262                DO i=1-OLx,sNx+OLx
263                 fVerP  (i,j,1,ip) = 0. _d 0
264  #ifdef ALLOW_AUTODIFF_TAMC               fVerP  (i,j,2,ip) = 0. _d 0
265  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte               kappaRTr(i,j,ip)  = 0. _d 0
266  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte              ENDDO
267  #ifdef ALLOW_KPP             ENDDO
268  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte           ENDDO
269  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  C-      set tracer tendency to zero:
270             DO iTracer=1,PTRACERS_num
271              DO k=1,Nr
272               DO j=1-OLy,sNy+OLy
273                DO i=1-OLx,sNx+OLx
274                 gPTr(i,j,k,bi,bj,itracer) = 0. _d 0
275                ENDDO
276               ENDDO
277              ENDDO
278             ENDDO
279            ENDIF
280  #endif  #endif
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 C--     Start of diagnostic loop  
         DO k=Nr,1,-1  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 C? Patrick, is this formula correct now that we change the loop range?  
 C? Do we still need this?  
 cph kkey formula corrected.  
 cph Needed for rhok, rhokm1, in the case useGMREDI.  
          kkey = (itdkey-1)*Nr + k  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 C--       Integrate continuity vertically for vertical velocity  
 c         CALL INTEGRATE_FOR_W(  
 c    I                         bi, bj, k, uVel, vVel,  
 c    O                         wVel,  
 c    I                         myThid )  
   
 #ifdef    ALLOW_OBCS  
 #ifdef    ALLOW_NONHYDROSTATIC  
 C--       Apply OBC to W if in N-H mode  
 c         IF (useOBCS.AND.nonHydrostatic) THEN  
 c           CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )  
 c         ENDIF  
 #endif    /* ALLOW_NONHYDROSTATIC */  
 #endif    /* ALLOW_OBCS */  
   
 C--       Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h  
 C--       MOST of THERMODYNAMICS will be disabled  
 #ifndef SINGLE_LAYER_MODE  
   
 C--       Calculate gradients of potential density for isoneutral  
 C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)  
 c         IF ( k.GT.1 .AND. (useGMRedi.OR.ivdc_kappa.NE.0.) ) THEN  
           IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.) ) THEN  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  
 CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  
 CADJ STORE pressure(:,:,k,bi,bj) =  
 CADJ &     comlev1_bibj_k, key=kkey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
             CALL FIND_RHO(  
      I        bi, bj, iMin, iMax, jMin, jMax, k, k,  
      I        theta, salt,  
      O        rhoK,  
      I        myThid )  
   
             IF (k.GT.1) THEN  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  
 CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  
 CADJ STORE pressure(:,:,k-1,bi,bj) =  
 CADJ &     comlev1_bibj_k, key=kkey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
              CALL FIND_RHO(  
      I        bi, bj, iMin, iMax, jMin, jMax, k-1, k,  
      I        theta, salt,  
      O        rhoKm1,  
      I        myThid )  
             ENDIF  
             CALL GRAD_SIGMA(  
      I             bi, bj, iMin, iMax, jMin, jMax, k,  
      I             rhoK, rhoKm1, rhoK,  
      O             sigmaX, sigmaY, sigmaR,  
      I             myThid )  
           ENDIF  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE rhok   (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte  
 CADJ STORE rhokm1 (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
 C--       Implicit Vertical Diffusion for Convection  
 c ==> should use sigmaR !!!  
           IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN  
             CALL CALC_IVDC(  
      I        bi, bj, iMin, iMax, jMin, jMax, k,  
      I        rhoKm1, rhoK,  
      U        ConvectCount, KappaRT, KappaRS,  
      I        myTime, myIter, myThid)  
           ENDIF  
   
 #endif /* SINGLE_LAYER_MODE */  
281    
282  C--     end of diagnostic k loop (Nr:1)  #ifdef ALLOW_ADAMSBASHFORTH_3
283          ENDDO  C-      Apply AB on T,S :
284            iterNb = myIter
285            IF (staggerTimeStep) iterNb = myIter - 1
286            m1 = 1 + MOD(iterNb+1,2)
287            m2 = 1 + MOD( iterNb ,2)
288    C       compute T^n+1/2 (stored in gtNm) extrapolating T forward in time
289            IF ( AdamsBashforth_T ) CALL ADAMS_BASHFORTH3(
290         I                                  bi, bj, 0,
291         U                                  theta, gtNm,
292         I                                  tempStartAB, iterNb, myThid )
293    C       compute S^n+1/2 (stored in gsNm) extrapolating S forward in time
294            IF ( AdamsBashforth_S ) CALL ADAMS_BASHFORTH3(
295         I                                  bi, bj, 0,
296         U                                  salt, gsNm,
297         I                                  saltStartAB, iterNb, myThid )
298    #endif /* ALLOW_ADAMSBASHFORTH_3 */
299    
300    c       iMin = 1-OLx
301    c       iMax = sNx+OLx
302    c       jMin = 1-OLy
303    c       jMax = sNy+OLy
304    
305  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
306  cph avoids recomputation of integrate_for_w  cph avoids recomputation of integrate_for_w
307  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
 CADJ STORE pressure (:,:,:,bi,bj) =  
 CADJ &     comlev1_bibj, key=itdkey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 #ifdef  ALLOW_OBCS  
 C--     Calculate future values on open boundaries  
         IF (useOBCS) THEN  
           CALL OBCS_CALC( bi, bj, myTime+deltaT, myIter+1,  
      I            uVel, vVel, wVel, theta, salt,  
      I            myThid )  
         ENDIF  
 #endif  /* ALLOW_OBCS */  
   
   
 c********************************************  
 cswdice --- add ---  
 #ifdef ALLOW_TSEAICE  
 C--     Determines forcing terms based on external fields  
 c--     including effects from ice  
         CALL ICE_FORCING(  
      I             bi, bj, iMin, iMax, jMin, jMax,  
      I             myThid )  
 #else  
   
 cswdice --- end add ---  
   
 C--     Determines forcing terms based on external fields  
 C       relaxation terms, etc.  
         CALL EXTERNAL_FORCING_SURF(  
      I             bi, bj, iMin, iMax, jMin, jMax,  
      I             myThid )  
 cswdice --- add ----  
 #endif  
 cswdice --- end add ---  
 c******************************************  
   
   
   
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 cph needed for KPP  
 CADJ STORE surfacetendencyU(:,:,bi,bj)  
 CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  
 CADJ STORE surfacetendencyV(:,:,bi,bj)  
 CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  
 CADJ STORE surfacetendencyS(:,:,bi,bj)  
 CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  
 CADJ STORE surfacetendencyT(:,:,bi,bj)  
 CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte  
308  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
309    
310  C--     Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h  C--     Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h
311  C--     MOST of THERMODYNAMICS will be disabled  C--     MOST of THERMODYNAMICS will be disabled
312  #ifndef SINGLE_LAYER_MODE  #ifndef SINGLE_LAYER_MODE
313    
 #ifdef  ALLOW_GMREDI  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE sigmaX(:,:,k)        = comlev1_bibj_k, key=kkey, byte=isbyte  
 CADJ STORE sigmaY(:,:,k)        = comlev1_bibj_k, key=kkey, byte=isbyte  
 CADJ STORE sigmaR(:,:,k)        = comlev1_bibj_k, key=kkey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
 C--     Calculate iso-neutral slopes for the GM/Redi parameterisation  
         IF (useGMRedi) THEN  
           CALL GMREDI_CALC_TENSOR(  
      I             bi, bj, iMin, iMax, jMin, jMax,  
      I             sigmaX, sigmaY, sigmaR,  
      I             myThid )  
 #ifdef ALLOW_AUTODIFF_TAMC  
         ELSE  
           CALL GMREDI_CALC_TENSOR_DUMMY(  
      I             bi, bj, iMin, iMax, jMin, jMax,  
      I             sigmaX, sigmaY, sigmaR,  
      I             myThid )  
 #endif /* ALLOW_AUTODIFF_TAMC */  
         ENDIF  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE Kwx(:,:,:,bi,bj)   = comlev1_bibj, key=itdkey, byte=isbyte  
 CADJ STORE Kwy(:,:,:,bi,bj)   = comlev1_bibj, key=itdkey, byte=isbyte  
 CADJ STORE Kwz(:,:,:,bi,bj)   = comlev1_bibj, key=itdkey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 #endif  /* ALLOW_GMREDI */  
   
 #ifdef  ALLOW_KPP  
 C--     Compute KPP mixing coefficients  
         IF (useKPP) THEN  
           CALL KPP_CALC(  
      I                  bi, bj, myTime, myThid )  
 #ifdef ALLOW_AUTODIFF_TAMC  
         ELSE  
           CALL KPP_CALC_DUMMY(  
      I                  bi, bj, myTime, myThid )  
 #endif /* ALLOW_AUTODIFF_TAMC */  
         ENDIF  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE KPPghat   (:,:,:,bi,bj)  
 CADJ &   , KPPdiffKzT(:,:,:,bi,bj)  
 CADJ &   , KPPdiffKzS(:,:,:,bi,bj)  
 CADJ &   , KPPfrac   (:,:  ,bi,bj)  
 CADJ &                 = comlev1_bibj, key=itdkey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 #endif  /* ALLOW_KPP */  
   
314  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
 CADJ STORE KappaRT(:,:,:)     = comlev1_bibj, key=itdkey, byte=isbyte  
 CADJ STORE KappaRS(:,:,:)     = comlev1_bibj, key=itdkey, byte=isbyte  
315  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
316  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
317  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
318  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
319  #ifdef ALLOW_PASSIVE_TRACER  cph-test
320  CADJ STORE tr1  (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  # ifdef ALLOW_DEPTH_CONTROL
321  #endif  CADJ STORE gtnm1(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
322    CADJ STORE gsnm1(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
323    # endif
324  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
325    
 #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( 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 */  
   
326  #ifndef DISABLE_MULTIDIM_ADVECTION  #ifndef DISABLE_MULTIDIM_ADVECTION
327  C--     Some advection schemes are better calculated using a multi-dimensional  C--     Some advection schemes are better calculated using a multi-dimensional
328  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 527  C recomputation. It *is* differentiable, Line 335  C recomputation. It *is* differentiable,
335  C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to  C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to
336  C disable this section of code.  C disable this section of code.
337          IF (tempMultiDimAdvec) THEN          IF (tempMultiDimAdvec) THEN
338            CALL GAD_ADVECTION(bi,bj,tempAdvScheme,GAD_TEMPERATURE,  #ifdef ALLOW_DEBUG
339       U                      theta,gT,            IF ( debugLevel .GE. debLevB )
340       I                      myTime,myIter,myThid)       &     CALL DEBUG_CALL('GAD_ADVECTION',myThid)
341    #endif
342              CALL GAD_ADVECTION(
343         I             tempImplVertAdv, tempAdvScheme, tempVertAdvScheme,
344         I             GAD_TEMPERATURE,
345         I             uVel, vVel, wVel, theta,
346         O             gT,
347         I             bi,bj,myTime,myIter,myThid)
348          ENDIF          ENDIF
349          IF (saltMultiDimAdvec) THEN          IF (saltMultiDimAdvec) THEN
350            CALL GAD_ADVECTION(bi,bj,saltAdvScheme,GAD_SALINITY,  #ifdef ALLOW_DEBUG
351       U                      salt,gS,            IF ( debugLevel .GE. debLevB )
352       I                      myTime,myIter,myThid)       &     CALL DEBUG_CALL('GAD_ADVECTION',myThid)
353    #endif
354              CALL GAD_ADVECTION(
355         I             saltImplVertAdv, saltAdvScheme, saltVertAdvScheme,
356         I             GAD_SALINITY,
357         I             uVel, vVel, wVel, salt,
358         O             gS,
359         I             bi,bj,myTime,myIter,myThid)
360          ENDIF          ENDIF
361    
362  C Since passive tracers are configurable separately from T,S we  C Since passive tracers are configurable separately from T,S we
363  C call the multi-dimensional method for PTRACERS regardless  C call the multi-dimensional method for PTRACERS regardless
364  C of whether multiDimAdvection is set or not.  C of whether multiDimAdvection is set or not.
365  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
366          IF ( usePTRACERS ) THEN          IF ( usePTRACERS ) THEN
367    #ifdef ALLOW_DEBUG
368              IF ( debugLevel .GE. debLevB )
369         &     CALL DEBUG_CALL('PTRACERS_ADVECTION',myThid)
370    #endif
371           CALL PTRACERS_ADVECTION( bi,bj,myIter,myTime,myThid )           CALL PTRACERS_ADVECTION( bi,bj,myIter,myTime,myThid )
372          ENDIF          ENDIF
373  #endif /* ALLOW_PTRACERS */  #endif /* ALLOW_PTRACERS */
374  #endif /* DISABLE_MULTIDIM_ADVECTION */  #endif /* DISABLE_MULTIDIM_ADVECTION */
375    
376    #ifdef ALLOW_DEBUG
377            IF ( debugLevel .GE. debLevB )
378         &    CALL DEBUG_MSG('ENTERING DOWNWARD K LOOP',myThid)
379    #endif
380    
381  C--     Start of thermodynamics loop  C--     Start of thermodynamics loop
382          DO k=Nr,1,-1          DO k=Nr,1,-1
383  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
384  C? Patrick Is this formula correct?  C? Patrick Is this formula correct?
385  cph Yes, but I rewrote it.  cph Yes, but I rewrote it.
386  cph Also, the KappaR? need the index and subscript k!  cph Also, the kappaR? need the index and subscript k!
387           kkey = (itdkey-1)*Nr + k           kkey = (itdkey-1)*Nr + k
388  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
389    
# Line 568  C--       kDown  Cycles through 2,1 to p Line 400  C--       kDown  Cycles through 2,1 to p
400            jMin = 1-OLy            jMin = 1-OLy
401            jMax = sNy+OLy            jMax = sNy+OLy
402    
403  C--       Get temporary terms used by tendency routines            IF (k.EQ.Nr) THEN
404               DO j=1-Oly,sNy+Oly
405                DO i=1-Olx,sNx+Olx
406                 rTransKp1(i,j) = 0. _d 0
407                ENDDO
408               ENDDO
409              ELSE
410               DO j=1-Oly,sNy+Oly
411                DO i=1-Olx,sNx+Olx
412                 rTransKp1(i,j) = rTrans(i,j)
413                ENDDO
414               ENDDO
415              ENDIF
416    #ifdef ALLOW_AUTODIFF_TAMC
417    CADJ STORE rTransKp1(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
418    #endif
419    
420    C--       Get temporary terms used by tendency routines :
421    C-        Calculate horizontal "volume transport" through tracer cell face
422            CALL CALC_COMMON_FACTORS (            CALL CALC_COMMON_FACTORS (
423       I         bi,bj,iMin,iMax,jMin,jMax,k,       I         uVel, vVel,
424       O         xA,yA,uTrans,vTrans,rTrans,maskUp,       O         uFld, vFld, uTrans, vTrans, xA, yA,
425       I         myThid)       I         k,bi,bj, myThid )
426    
427    C-        Calculate vertical "volume transport" through tracer cell face
428              IF (k.EQ.1) THEN
429    C-         Surface interface :
430               DO j=1-Oly,sNy+Oly
431                DO i=1-Olx,sNx+Olx
432                 wFld(i,j)   = 0. _d 0
433                 maskUp(i,j) = 0. _d 0
434                 rTrans(i,j) = 0. _d 0
435                ENDDO
436               ENDDO
437              ELSE
438    C-         Interior interface :
439               DO j=1-Oly,sNy+Oly
440                DO i=1-Olx,sNx+Olx
441                 wFld(i,j)   = wVel(i,j,k,bi,bj)
442                 maskUp(i,j) = maskC(i,j,k-1,bi,bj)*maskC(i,j,k,bi,bj)
443                 rTrans(i,j) = wFld(i,j)*rA(i,j,bi,bj)*maskUp(i,j)
444                ENDDO
445               ENDDO
446              ENDIF
447    
448  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
449  C--   Residual transp = Bolus transp + Eulerian transp  C--   Residual transp = Bolus transp + Eulerian transp
450            IF (useGMRedi) THEN            IF (useGMRedi) THEN
451              CALL GMREDI_CALC_UVFLOW(              CALL GMREDI_CALC_UVFLOW(
452       &            uTrans, vTrans, bi, bj, k, myThid)       U                  uFld, vFld, uTrans, vTrans,
453         I                  k, bi, bj, myThid )
454              IF (K.GE.2) CALL GMREDI_CALC_WFLOW(              IF (K.GE.2) CALL GMREDI_CALC_WFLOW(
455       &                    rTrans, bi, bj, k, myThid)       U                  wFld, rTrans,
456         I                  k, bi, bj, myThid )
457            ENDIF            ENDIF
458    # ifdef ALLOW_AUTODIFF_TAMC
459    CADJ STORE rTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
460    # ifdef GM_BOLUS_ADVEC
461    CADJ STORE uTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
462    CADJ STORE vTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
463    # endif
464    # endif /* ALLOW_AUTODIFF_TAMC */
465  #endif /* ALLOW_GMREDI */  #endif /* ALLOW_GMREDI */
466    
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE KappaRT(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte  
 CADJ STORE KappaRS(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
467  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL
468  C--      Calculate the total vertical diffusivity  C--      Calculate the total vertical diffusivity
469           CALL CALC_DIFFUSIVITY(            IF ( .NOT.implicitDiffusion ) THEN
470       I        bi,bj,iMin,iMax,jMin,jMax,k,              CALL CALC_DIFFUSIVITY(
471       I        maskUp,       I          bi,bj,iMin,iMax,jMin,jMax,k,
472       O        KappaRT,KappaRS,       I          maskUp,
473       I        myThid)       O          kappaRT,kappaRS,
474         I          myThid)
475              ENDIF
476    # ifdef ALLOW_AUTODIFF_TAMC
477    CADJ STORE kappaRT(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
478    CADJ STORE kappaRS(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
479    # endif /* ALLOW_AUTODIFF_TAMC */
480  #endif  #endif
481    
482            iMin = 1-OLx+2            iMin = 1-OLx+2
# Line 604  C--      Calculate the total vertical di Line 485  C--      Calculate the total vertical di
485            jMax = sNy+OLy-1            jMax = sNy+OLy-1
486    
487  C--      Calculate active tracer tendencies (gT,gS,...)  C--      Calculate active tracer tendencies (gT,gS,...)
488  C        and step forward storing result in gTnm1, gSnm1, etc.  C        and step forward storing result in gT, gS, etc.
489    C--
490    # ifdef ALLOW_AUTODIFF_TAMC
491    #  if ((defined NONLIN_FRSURF) || (defined ALLOW_DEPTH_CONTROL)) && (defined ALLOW_GMREDI)
492    CADJ STORE kux(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
493    CADJ STORE kvy(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
494    #   ifdef GM_EXTRA_DIAGONAL
495    CADJ STORE kuz(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
496    CADJ STORE kvz(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
497    #   endif
498    #  endif
499    # endif /* ALLOW_AUTODIFF_TAMC */
500           IF ( tempStepping ) THEN           IF ( tempStepping ) THEN
501             CALL CALC_GT(  #ifdef ALLOW_AUTODIFF_TAMC
502    CADJ STORE gTnm1(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
503    # if (defined NONLIN_FRSURF) || (defined ALLOW_DEPTH_CONTROL)
504    CADJ STORE gT(:,:,k,bi,bj)    = comlev1_bibj_k, key=kkey, byte=isbyte
505    # endif
506    #endif /* ALLOW_AUTODIFF_TAMC */
507              CALL CALC_GT(
508       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
509       I         xA,yA,uTrans,vTrans,rTrans,maskUp,       I         xA, yA, maskUp, uFld, vFld, wFld,
510       I         KappaRT,       I         uTrans, vTrans, rTrans, rTransKp1,
511         I         kappaRT,
512       U         fVerT,       U         fVerT,
513       I         myTime,myIter,myThid)       I         myTime,myIter,myThid)
514    #ifdef ALLOW_ADAMSBASHFORTH_3
515              IF ( AdamsBashforth_T ) THEN
516               CALL TIMESTEP_TRACER(
517         I         bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,
518         I         gtNm(1-Olx,1-Oly,1,1,1,m2),
519         U         gT,
520         I         myIter, myThid)
521              ELSE
522    #endif
523             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
524       I         bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,       I         bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,
525       I         theta, gT,       I         theta,
526         U         gT,
527       I         myIter, myThid)       I         myIter, myThid)
528    #ifdef ALLOW_ADAMSBASHFORTH_3
529              ENDIF
530    #endif
531           ENDIF           ENDIF
532    
533           IF ( saltStepping ) THEN           IF ( saltStepping ) THEN
534             CALL CALC_GS(  #ifdef ALLOW_AUTODIFF_TAMC
535    CADJ STORE gSnm1(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
536    # if (defined NONLIN_FRSURF) || (defined ALLOW_DEPTH_CONTROL)
537    CADJ STORE gS(:,:,k,bi,bj)    = comlev1_bibj_k, key=kkey, byte=isbyte
538    # endif
539    #endif /* ALLOW_AUTODIFF_TAMC */
540              CALL CALC_GS(
541       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
542       I         xA,yA,uTrans,vTrans,rTrans,maskUp,       I         xA, yA, maskUp, uFld, vFld, wFld,
543       I         KappaRS,       I         uTrans, vTrans, rTrans, rTransKp1,
544         I         kappaRS,
545       U         fVerS,       U         fVerS,
546       I         myTime,myIter,myThid)       I         myTime,myIter,myThid)
547    #ifdef ALLOW_ADAMSBASHFORTH_3
548              IF ( AdamsBashforth_S ) THEN
549             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
550       I         bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,       I         bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,
551       I         salt, gS,       I         gsNm(1-Olx,1-Oly,1,1,1,m2),
552         U         gS,
553       I         myIter, myThid)       I         myIter, myThid)
554           ENDIF            ELSE
555  #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)  
556             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
557       I         bi,bj,iMin,iMax,jMin,jMax,k,tracerAdvScheme,       I         bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,
558       I         Tr1, gTr1,       I         salt,
559       I         myIter,myThid)       U         gS,
560           ENDIF       I         myIter, myThid)
561    #ifdef ALLOW_ADAMSBASHFORTH_3
562              ENDIF
563  #endif  #endif
564             ENDIF
565    
566  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
567           IF ( usePTRACERS ) THEN           IF ( usePTRACERS ) THEN
568             CALL PTRACERS_INTEGERATE(             IF ( .NOT.implicitDiffusion ) THEN
569       I         bi,bj,k,               CALL PTRACERS_CALC_DIFF(
570       I         xA,yA,uTrans,vTrans,rTrans,maskUp,       I            bi,bj,iMin,iMax,jMin,jMax,k,
571       X         KappaRS,       I            maskUp,
572       I         myIter,myTime,myThid)       O            kappaRTr,
573         I            myThid)
574               ENDIF
575    # ifdef ALLOW_AUTODIFF_TAMC
576    CADJ STORE kappaRTr(:,:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
577    # endif /* ALLOW_AUTODIFF_TAMC */
578               CALL PTRACERS_INTEGRATE(
579         I          bi,bj,k,
580         I          xA, yA, maskUp, uFld, vFld, wFld,
581         I          uTrans, vTrans, rTrans, rTransKp1,
582         I          kappaRTr,
583         U          fVerP,
584         I          myTime,myIter,myThid)
585           ENDIF           ENDIF
586  #endif /* ALLOW_PTRACERS */  #endif /* ALLOW_PTRACERS */
587    
# Line 661  C--      Apply open boundary conditions Line 593  C--      Apply open boundary conditions
593  #endif   /* ALLOW_OBCS */  #endif   /* ALLOW_OBCS */
594    
595  C--      Freeze water  C--      Freeze water
596           IF ( allowFreezing .AND. .NOT. useSEAICE ) THEN  C  this bit of code is left here for backward compatibility.
597    C  freezing at surface level has been moved to FORWARD_STEP
598             IF ( useOldFreezing .AND. .NOT. useSEAICE
599         &       .AND. .NOT.(useThSIce.AND.k.EQ.1) ) THEN
600  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
601  CADJ STORE gT(:,:,k,bi,bj) = comlev1_bibj_k  CADJ STORE gT(:,:,k,bi,bj) = comlev1_bibj_k
602  CADJ &   , key = kkey, byte = isbyte  CADJ &   , key = kkey, byte = isbyte
603  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
604              CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )              CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )
605           END IF           ENDIF
606    
607  C--     end of thermodynamic k loop (Nr:1)  C--     end of thermodynamic k loop (Nr:1)
608          ENDDO          ENDDO
609    
610  cswdice -- add ---  C       All explicit advection/diffusion/sources should now be
611  #ifdef ALLOW_TSEAICE  C       done. The updated tracer field is in gPtr. Accumalate
612  c timeaveraging for ice model values  C       explicit tendency and also reset gPtr to initial tracer
613             CALL ICE_AVE(bi,bj,iMin,iMax,jMin,jMax,myThid )  C       field for implicit matrix calculation
614  #endif  
615  cswdice --- end add ---  #ifdef ALLOW_MATRIX
616            IF (useMATRIX)
617         &    CALL MATRIX_STORE_TENDENCY_EXP(bi,bj, myTime,myIter,myThid)
618    #endif
619    
620  C--     Implicit diffusion          iMin = 1
621          IF (implicitDiffusion) THEN          iMax = sNx
622            jMin = 1
623           IF (tempStepping) THEN          jMax = sNy
624    
625    C--     Implicit vertical advection & diffusion
626            IF ( tempStepping .AND. implicitDiffusion ) THEN
627              CALL CALC_3D_DIFFUSIVITY(
628         I         bi,bj,iMin,iMax,jMin,jMax,
629         I         GAD_TEMPERATURE, useGMredi, useKPP,
630         O         kappaRk,
631         I         myThid)
632            ENDIF
633    #ifdef INCLUDE_IMPLVERTADV_CODE
634            IF ( tempImplVertAdv ) THEN
635              CALL GAD_IMPLICIT_R(
636         I         tempImplVertAdv, tempAdvScheme, GAD_TEMPERATURE,
637         I         kappaRk, wVel, theta,
638         U         gT,
639         I         bi, bj, myTime, myIter, myThid )
640            ELSEIF ( tempStepping .AND. implicitDiffusion ) THEN
641    #else /* INCLUDE_IMPLVERTADV_CODE */
642            IF     ( tempStepping .AND. implicitDiffusion ) THEN
643    #endif /* INCLUDE_IMPLVERTADV_CODE */
644  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
645    CADJ STORE kappaRk(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
646  CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte  CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
647  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
648              CALL IMPLDIFF(            CALL IMPLDIFF(
649       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
650       I         deltaTtracer, KappaRT, recip_HFacC,       I         GAD_TEMPERATURE, kappaRk, recip_hFacC,
651       U         gT,       U         gT,
652       I         myThid )       I         myThid )
653            ENDIF
654    
655    #ifdef ALLOW_TIMEAVE
656            useVariableK = useKPP .OR. usePP81 .OR. useMY82 .OR. useGGL90
657         &       .OR. useGMredi .OR. ivdc_kappa.NE.0.
658            IF (taveFreq.GT.0. .AND. useVariableK ) THEN
659             IF (implicitDiffusion) THEN
660               CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, gT, kappaRk,
661         I                        Nr, 3, deltaTclock, bi, bj, myThid)
662    c        ELSE
663    c          CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, theta, kappaRT,
664    c    I                        Nr, 3, deltaTclock, bi, bj, myThid)
665           ENDIF           ENDIF
666            ENDIF
667    #endif /* ALLOW_TIMEAVE */
668    
669            IF ( saltStepping .AND. implicitDiffusion ) THEN
670              CALL CALC_3D_DIFFUSIVITY(
671         I         bi,bj,iMin,iMax,jMin,jMax,
672         I         GAD_SALINITY, useGMredi, useKPP,
673         O         kappaRk,
674         I         myThid)
675            ENDIF
676    
677           IF (saltStepping) THEN  #ifdef INCLUDE_IMPLVERTADV_CODE
678            IF ( saltImplVertAdv ) THEN
679              CALL GAD_IMPLICIT_R(
680         I         saltImplVertAdv, saltAdvScheme, GAD_SALINITY,
681         I         kappaRk, wVel, salt,
682         U         gS,
683         I         bi, bj, myTime, myIter, myThid )
684            ELSEIF ( saltStepping .AND. implicitDiffusion ) THEN
685    #else /* INCLUDE_IMPLVERTADV_CODE */
686            IF     ( saltStepping .AND. implicitDiffusion ) THEN
687    #endif /* INCLUDE_IMPLVERTADV_CODE */
688  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
689    CADJ STORE kappaRk(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
690  CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte  CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
691  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
692              CALL IMPLDIFF(            CALL IMPLDIFF(
693       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
694       I         deltaTtracer, KappaRS, recip_HFacC,       I         GAD_SALINITY, kappaRk, recip_hFacC,
695       U         gS,       U         gS,
696       I         myThid )       I         myThid )
697           ENDIF          ENDIF
   
 #ifdef ALLOW_PASSIVE_TRACER  
          IF (tr1Stepping) THEN  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE gTr1(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, 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  
698    
699  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
700  C Vertical diffusion (implicit) for passive tracers          IF     ( usePTRACERS ) THEN
701           IF ( usePTRACERS ) THEN  C--     Vertical advection/diffusion (implicit) for passive tracers
702             CALL PTRACERS_IMPLDIFF( bi,bj,KappaRS,myThid )             CALL PTRACERS_IMPLICIT(
703           ENDIF       U                             kappaRk,
704         I                             bi, bj, myTime, myIter, myThid )
705            ENDIF
706  #endif /* ALLOW_PTRACERS */  #endif /* ALLOW_PTRACERS */
707    
708  #ifdef   ALLOW_OBCS  #ifdef   ALLOW_OBCS
709  C--      Apply open boundary conditions  C--      Apply open boundary conditions
710           IF (useOBCS) THEN          IF ( ( implicitDiffusion
711         &    .OR. tempImplVertAdv
712         &    .OR. saltImplVertAdv
713         &       ) .AND. useOBCS     ) THEN
714             DO K=1,Nr             DO K=1,Nr
715               CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )               CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )
716             ENDDO             ENDDO
          END IF  
 #endif   /* ALLOW_OBCS */  
   
 C--     End If implicitDiffusion  
717          ENDIF          ENDIF
718    #endif   /* ALLOW_OBCS */
719    
720  #endif /* SINGLE_LAYER_MODE */  #endif /* SINGLE_LAYER_MODE */
721    
722  Ccs-  C--   end bi,bj loops.
723         ENDDO         ENDDO
724        ENDDO        ENDDO
725    
726  #ifdef ALLOW_AIM  #ifdef ALLOW_DEBUG
 c     IF ( useAIM ) THEN  
 c      CALL AIM_AIM2DYN_EXCHANGES( myTime, myThid )  
 c     ENDIF  
 #endif /* ALLOW_AIM */  
 c     IF ( staggerTimeStep ) THEN  
 c      IF ( useAIM .OR. useCubedSphereExchange ) THEN  
 c        IF (tempStepping) _EXCH_XYZ_R8(gT,myThid)  
 c        IF (saltStepping) _EXCH_XYZ_R8(gS,myThid)  
 c      ELSEIF ( useGMRedi .AND. Oly.LT.4 ) THEN  
 cc       .AND. GM_AdvForm .AND. .NOT.GM_AdvSeparate ) THEN  
 c        IF (tempMultiDimAdvec) _EXCH_XYZ_R8(gT,myThid)  
 c        IF (saltMultiDimAdvec) _EXCH_XYZ_R8(gS,myThid)  
 c      ENDIF    
 c     ENDIF    
   
 #ifndef DISABLE_DEBUGMODE  
727        If (debugMode) THEN        If (debugMode) THEN
728         CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (THERMODYNAMICS)',myThid)
729         CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (THERMODYNAMICS)',myThid)
730         CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (THERMODYNAMICS)',myThid)
731         CALL DEBUG_STATS_RL(Nr,theta,'Theta (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,theta,'Theta (THERMODYNAMICS)',myThid)
732         CALL DEBUG_STATS_RL(Nr,salt,'Salt (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,salt,'Salt (THERMODYNAMICS)',myThid)
733         CALL DEBUG_STATS_RL(Nr,Gt,'Gt (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,gT,'Gt (THERMODYNAMICS)',myThid)
734         CALL DEBUG_STATS_RL(Nr,Gs,'Gs (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,gS,'Gs (THERMODYNAMICS)',myThid)
735         CALL DEBUG_STATS_RL(Nr,GtNm1,'GtNm1 (THERMODYNAMICS)',myThid)  #ifndef ALLOW_ADAMSBASHFORTH_3
736         CALL DEBUG_STATS_RL(Nr,GsNm1,'GsNm1 (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,gtNm1,'GtNm1 (THERMODYNAMICS)',myThid)
737           CALL DEBUG_STATS_RL(Nr,gsNm1,'GsNm1 (THERMODYNAMICS)',myThid)
738    #endif
739  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
740         IF ( usePTRACERS ) THEN         IF ( usePTRACERS ) THEN
741           CALL PTRACERS_DEBUG(myThid)           CALL PTRACERS_DEBUG(myThid)
742         ENDIF         ENDIF
743  #endif /* ALLOW_PTRACERS */  #endif /* ALLOW_PTRACERS */
744        ENDIF        ENDIF
745    #endif /* ALLOW_DEBUG */
746    
747    #ifdef ALLOW_DEBUG
748             IF ( debugLevel .GE. debLevB )
749         &    CALL DEBUG_LEAVE('THERMODYNAMICS',myThid)
750  #endif  #endif
751    
752    #endif /* ALLOW_GENERIC_ADVDIFF */
753    
754        RETURN        RETURN
755        END        END
Legend:
Removed from v.1.31  
changed lines
  Added in v.1.105
  ViewVC Help
Powered by ViewVC 1.1.22