/[MITgcm]/MITgcm/model/src/thermodynamics.F
ViewVC logotype

Diff of /MITgcm/model/src/thermodynamics.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.24 by jmc, Wed Jul 3 20:22:39 2002 UTC revision 1.104 by heimbach, Wed Jun 7 01:55:13 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 17  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 75  C     == Global variables === Line 77  C     == Global variables ===
77  #include "PARAMS.h"  #include "PARAMS.h"
78  #include "DYNVARS.h"  #include "DYNVARS.h"
79  #include "GRID.h"  #include "GRID.h"
80    #ifdef ALLOW_GENERIC_ADVDIFF
81  #include "GAD.h"  #include "GAD.h"
 #ifdef ALLOW_PASSIVE_TRACER  
 #include "TR1.h"  
82  #endif  #endif
83    #ifdef ALLOW_PTRACERS
84    #include "PTRACERS_SIZE.h"
85    #include "PTRACERS.h"
86    #endif
87    #ifdef ALLOW_TIMEAVE
88    #include "TIMEAVE_STATV.h"
89    #endif
90    
91  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
92  # include "tamc.h"  # include "tamc.h"
93  # include "tamc_keys.h"  # include "tamc_keys.h"
94  # include "FFIELDS.h"  # include "FFIELDS.h"
95    # include "EOS.h"
96  # ifdef ALLOW_KPP  # ifdef ALLOW_KPP
97  #  include "KPP.h"  #  include "KPP.h"
98  # endif  # endif
99  # ifdef ALLOW_GMREDI  # ifdef ALLOW_GMREDI
100  #  include "GMREDI.h"  #  include "GMREDI.h"
101  # endif  # endif
102    # ifdef ALLOW_EBM
103    #  include "EBM.h"
104    # endif
105  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
 #ifdef ALLOW_TIMEAVE  
 #include "TIMEAVE_STATV.h"  
 #endif  
106    
107  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
108  C     == Routine arguments ==  C     == Routine arguments ==
# Line 103  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
# Line 111  C                              transport Line 122  C                              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 137  C                      index into fVerTe Line 145  C                      index into fVerTe
145        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
146        _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
147        _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
148          _RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
149        _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
150        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
151        _RL fVerS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
152        _RL fVerTr1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL kappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
153        _RL phiHyd  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL kappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
154        _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  #ifdef ALLOW_PTRACERS
155        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fVerP   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2,PTRACERS_num)
156        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL kappaRTr(1-Olx:sNx+Olx,1-Oly:sNy+Oly,PTRACERS_num)
157        _RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  #endif
158        _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)  
159        INTEGER iMin, iMax        INTEGER iMin, iMax
160        INTEGER jMin, jMax        INTEGER jMin, jMax
161        INTEGER bi, bj        INTEGER bi, bj
162        INTEGER i, j        INTEGER i, j
163        INTEGER k, km1, kup, kDown        INTEGER k, km1, kup, kDown
164    #ifdef ALLOW_ADAMSBASHFORTH_3
165          INTEGER iterNb, m1, m2
166    #endif
167    #ifdef ALLOW_TIMEAVE
168          LOGICAL useVariableK
169    #endif
170    #ifdef ALLOW_PTRACERS
171          INTEGER iTracer, ip
172    #endif
173    
174  CEOP  CEOP
175    
176    #ifdef ALLOW_DEBUG
177             IF ( debugLevel .GE. debLevB )
178         &    CALL DEBUG_ENTER('THERMODYNAMICS',myThid)
179    #endif
180    
181  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
182  C--   dummy statement to end declaration part  C--   dummy statement to end declaration part
183        ikey = 1        ikey = 1
184          itdkey = 1
185  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
186    
 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  
   
   
187  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
188  C--   HPF directive to help TAMC  C--   HPF directive to help TAMC
189  CHPF$ INDEPENDENT  CHPF$ INDEPENDENT
# Line 194  CHPF$ INDEPENDENT Line 194  CHPF$ INDEPENDENT
194  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
195  C--    HPF directive to help TAMC  C--    HPF directive to help TAMC
196  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS
197  CHPF$&                  ,phiHyd,utrans,vtrans,xA,yA  CHPF$&                  ,utrans,vtrans,xA,yA
198  CHPF$&                  ,KappaRT,KappaRS  CHPF$&                  ,kappaRT,kappaRS
199  CHPF$&                  )  CHPF$&                  )
200    # ifdef ALLOW_PTRACERS
201    CHPF$  INDEPENDENT, NEW (fVerP,kappaRTr)
202    # endif
203  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
204    
205         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
# Line 209  CHPF$&                  ) Line 212  CHPF$&                  )
212            act3 = myThid - 1            act3 = myThid - 1
213            max3 = nTx*nTy            max3 = nTx*nTy
214            act4 = ikey_dynamics - 1            act4 = ikey_dynamics - 1
215            ikey = (act1 + 1) + act2*max1            itdkey = (act1 + 1) + act2*max1
216       &                      + act3*max1*max2       &                      + act3*max1*max2
217       &                      + act4*max1*max2*max3       &                      + act4*max1*max2*max3
218  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
219    
220  C--     Set up work arrays that need valid initial values  C--   Set up work arrays with valid (i.e. not NaN) values
221    C     These inital values do not alter the numerical results. They
222    C     just ensure that all memory references are to valid floating
223    C     point numbers. This prevents spurious hardware signals due to
224    C     uninitialised but inert locations.
225    
226          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
227           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
228              xA(i,j)        = 0. _d 0
229              yA(i,j)        = 0. _d 0
230              uTrans(i,j)    = 0. _d 0
231              vTrans(i,j)    = 0. _d 0
232            rTrans (i,j)   = 0. _d 0            rTrans (i,j)   = 0. _d 0
233              rTransKp1(i,j) = 0. _d 0
234            fVerT  (i,j,1) = 0. _d 0            fVerT  (i,j,1) = 0. _d 0
235            fVerT  (i,j,2) = 0. _d 0            fVerT  (i,j,2) = 0. _d 0
236            fVerS  (i,j,1) = 0. _d 0            fVerS  (i,j,1) = 0. _d 0
237            fVerS  (i,j,2) = 0. _d 0            fVerS  (i,j,2) = 0. _d 0
238            fVerTr1(i,j,1) = 0. _d 0            kappaRT(i,j)   = 0. _d 0
239            fVerTr1(i,j,2) = 0. _d 0            kappaRS(i,j)   = 0. _d 0
           rhoKM1 (i,j)   = 0. _d 0  
240           ENDDO           ENDDO
241          ENDDO          ENDDO
242    
# Line 232  C--     Set up work arrays that need val Line 244  C--     Set up work arrays that need val
244           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
245            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
246  C This is currently also used by IVDC and Diagnostics  C This is currently also used by IVDC and Diagnostics
247             phiHyd(i,j,k)  = 0. _d 0             kappaRk(i,j,k)    = 0. _d 0
248             sigmaX(i,j,k) = 0. _d 0  C- tracer tendency needs to be set to zero (moved here from gad_calc_rhs):
249             sigmaY(i,j,k) = 0. _d 0             gT(i,j,k,bi,bj)   = 0. _d 0
250             sigmaR(i,j,k) = 0. _d 0             gS(i,j,k,bi,bj)   = 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  
            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  
 #ifdef ALLOW_GMREDI  
            Kwx(i,j,k,bi,bj)    = 0. _d 0  
            Kwy(i,j,k,bi,bj)    = 0. _d 0  
            Kwz(i,j,k,bi,bj)    = 0. _d 0  
 #ifdef GM_NON_UNITY_DIAGONAL  
            Kux(i,j,k,bi,bj)    = 0. _d 0  
            Kvy(i,j,k,bi,bj)    = 0. _d 0  
 #endif  
 #endif /* ALLOW_GMREDI */  
 #endif  
251            ENDDO            ENDDO
252           ENDDO           ENDDO
253          ENDDO          ENDDO
254    
255          iMin = 1-OLx  #ifdef ALLOW_PTRACERS
256          iMax = sNx+OLx          IF ( usePTRACERS ) THEN
257          jMin = 1-OLy           DO ip=1,PTRACERS_num
258          jMax = sNy+OLy             DO j=1-OLy,sNy+OLy
259                DO i=1-OLx,sNx+OLx
260                 fVerP  (i,j,1,ip) = 0. _d 0
261  #ifdef ALLOW_AUTODIFF_TAMC               fVerP  (i,j,2,ip) = 0. _d 0
262  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte               kappaRTr(i,j,ip)  = 0. _d 0
263  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte              ENDDO
264  #ifdef ALLOW_KPP             ENDDO
265  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte           ENDDO
266  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte  C-      set tracer tendency to zero:
267             DO iTracer=1,PTRACERS_num
268              DO k=1,Nr
269               DO j=1-OLy,sNy+OLy
270                DO i=1-OLx,sNx+OLx
271                 gPTr(i,j,k,bi,bj,itracer) = 0. _d 0
272                ENDDO
273               ENDDO
274              ENDDO
275             ENDDO
276            ENDIF
277  #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 = (ikey-1)*Nr + k  
 CADJ STORE rhokm1(:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte  
 CADJ STORE rhok  (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 C--       Integrate continuity vertically for vertical velocity  
           CALL INTEGRATE_FOR_W(  
      I                         bi, bj, k, uVel, vVel,  
      O                         wVel,  
      I                         myThid )  
   
 #ifdef    ALLOW_OBCS  
 #ifdef    ALLOW_NONHYDROSTATIC  
 C--       Apply OBC to W if in N-H mode  
           IF (useOBCS.AND.nonHydrostatic) THEN  
             CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )  
           ENDIF  
 #endif    /* ALLOW_NONHYDROSTATIC */  
 #endif    /* ALLOW_OBCS */  
278    
279  C--       Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h  #ifdef ALLOW_ADAMSBASHFORTH_3
280  C--       MOST of THERMODYNAMICS will be disabled  C-      Apply AB on T,S :
281  #ifndef SINGLE_LAYER_MODE          iterNb = myIter
282            IF (staggerTimeStep) iterNb = myIter - 1
283  C--       Calculate gradients of potential density for isoneutral          m1 = 1 + MOD(iterNb+1,2)
284  C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)          m2 = 1 + MOD( iterNb ,2)
285  c         IF ( k.GT.1 .AND. (useGMRedi.OR.ivdc_kappa.NE.0.) ) THEN  C       compute T^n+1/2 (stored in gtNm) extrapolating T forward in time
286            IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.) ) THEN          IF ( AdamsBashforth_T ) CALL ADAMS_BASHFORTH3(
287  #ifdef ALLOW_AUTODIFF_TAMC       I                                  bi, bj, 0,
288  CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte       U                                  theta, gtNm,
289  CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte       I                                  tempStartAB, iterNb, myThid )
290  #endif /* ALLOW_AUTODIFF_TAMC */  C       compute S^n+1/2 (stored in gsNm) extrapolating S forward in time
291              CALL FIND_RHO(          IF ( AdamsBashforth_S ) CALL ADAMS_BASHFORTH3(
292       I        bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,       I                                  bi, bj, 0,
293       I        theta, salt,       U                                  salt, gsNm,
294       O        rhoK,       I                                  saltStartAB, iterNb, myThid )
295       I        myThid )  #endif /* ALLOW_ADAMSBASHFORTH_3 */
296              IF (k.GT.1) THEN  
297  #ifdef ALLOW_AUTODIFF_TAMC  c       iMin = 1-OLx
298  CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  c       iMax = sNx+OLx
299  CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  c       jMin = 1-OLy
300  #endif /* ALLOW_AUTODIFF_TAMC */  c       jMax = sNy+OLy
              CALL FIND_RHO(  
      I        bi, bj, iMin, iMax, jMin, jMax, k-1, k, eosType,  
      I        theta, salt,  
      O        rhoKm1,  
      I        myThid )  
             ENDIF  
             CALL GRAD_SIGMA(  
      I             bi, bj, iMin, iMax, jMin, jMax, k,  
      I             rhoK, rhoKm1, rhoK,  
      O             sigmaX, sigmaY, sigmaR,  
      I             myThid )  
           ENDIF  
   
 C--       Implicit Vertical Diffusion for Convection  
 c ==> should use sigmaR !!!  
           IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN  
             CALL CALC_IVDC(  
      I        bi, bj, iMin, iMax, jMin, jMax, k,  
      I        rhoKm1, rhoK,  
      U        ConvectCount, KappaRT, KappaRS,  
      I        myTime, myIter, myThid)  
           ENDIF  
   
 #endif /* SINGLE_LAYER_MODE */  
   
 C--     end of diagnostic k loop (Nr:1)  
         ENDDO  
301    
302  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
303  cph avoids recomputation of integrate_for_w  cph avoids recomputation of integrate_for_w
304  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  
305  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
306    
307  C--     Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h  C--     Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h
308  C--     MOST of THERMODYNAMICS will be disabled  C--     MOST of THERMODYNAMICS will be disabled
309  #ifndef SINGLE_LAYER_MODE  #ifndef SINGLE_LAYER_MODE
310    
 #ifdef  ALLOW_GMREDI  
   
311  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
312  CADJ STORE sigmaX(:,:,k)        = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
313  CADJ STORE sigmaY(:,:,k)        = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
314  CADJ STORE sigmaR(:,:,k)        = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
315  #endif /* ALLOW_AUTODIFF_TAMC */  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
316  C--     Calculate iso-neutral slopes for the GM/Redi parameterisation  cph-test
317          IF (useGMRedi) THEN  # ifdef ALLOW_DEPTH_CONTROL
318            CALL GMREDI_CALC_TENSOR(  CADJ STORE gtnm1(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
319       I             bi, bj, iMin, iMax, jMin, jMax,  CADJ STORE gsnm1(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
320       I             sigmaX, sigmaY, sigmaR,  # endif
      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  
321  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
322    
 #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 */  
   
323  #ifndef DISABLE_MULTIDIM_ADVECTION  #ifndef DISABLE_MULTIDIM_ADVECTION
324  C--     Some advection schemes are better calculated using a multi-dimensional  C--     Some advection schemes are better calculated using a multi-dimensional
325  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 480  C recomputation. It *is* differentiable, Line 332  C recomputation. It *is* differentiable,
332  C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to  C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to
333  C disable this section of code.  C disable this section of code.
334          IF (tempMultiDimAdvec) THEN          IF (tempMultiDimAdvec) THEN
335            CALL GAD_ADVECTION(bi,bj,tempAdvScheme,GAD_TEMPERATURE,  #ifdef ALLOW_DEBUG
336       U                      theta,gT,            IF ( debugLevel .GE. debLevB )
337       I                      myTime,myIter,myThid)       &     CALL DEBUG_CALL('GAD_ADVECTION',myThid)
338    #endif
339              CALL GAD_ADVECTION(
340         I             tempImplVertAdv, tempAdvScheme, tempVertAdvScheme,
341         I             GAD_TEMPERATURE,
342         I             uVel, vVel, wVel, theta,
343         O             gT,
344         I             bi,bj,myTime,myIter,myThid)
345          ENDIF          ENDIF
346          IF (saltMultiDimAdvec) THEN          IF (saltMultiDimAdvec) THEN
347            CALL GAD_ADVECTION(bi,bj,saltAdvScheme,GAD_SALINITY,  #ifdef ALLOW_DEBUG
348       U                      salt,gS,            IF ( debugLevel .GE. debLevB )
349       I                      myTime,myIter,myThid)       &     CALL DEBUG_CALL('GAD_ADVECTION',myThid)
350    #endif
351              CALL GAD_ADVECTION(
352         I             saltImplVertAdv, saltAdvScheme, saltVertAdvScheme,
353         I             GAD_SALINITY,
354         I             uVel, vVel, wVel, salt,
355         O             gS,
356         I             bi,bj,myTime,myIter,myThid)
357          ENDIF          ENDIF
358    
359  C Since passive tracers are configurable separately from T,S we  C Since passive tracers are configurable separately from T,S we
360  C call the multi-dimensional method for PTRACERS regardless  C call the multi-dimensional method for PTRACERS regardless
361  C of whether multiDimAdvection is set or not.  C of whether multiDimAdvection is set or not.
362  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
363          IF ( usePTRACERS ) THEN          IF ( usePTRACERS ) THEN
364    #ifdef ALLOW_DEBUG
365              IF ( debugLevel .GE. debLevB )
366         &     CALL DEBUG_CALL('PTRACERS_ADVECTION',myThid)
367    #endif
368           CALL PTRACERS_ADVECTION( bi,bj,myIter,myTime,myThid )           CALL PTRACERS_ADVECTION( bi,bj,myIter,myTime,myThid )
369          ENDIF          ENDIF
370  #endif /* ALLOW_PTRACERS */  #endif /* ALLOW_PTRACERS */
371  #endif /* DISABLE_MULTIDIM_ADVECTION */  #endif /* DISABLE_MULTIDIM_ADVECTION */
372    
373    #ifdef ALLOW_DEBUG
374            IF ( debugLevel .GE. debLevB )
375         &    CALL DEBUG_MSG('ENTERING DOWNWARD K LOOP',myThid)
376    #endif
377    
378  C--     Start of thermodynamics loop  C--     Start of thermodynamics loop
379          DO k=Nr,1,-1          DO k=Nr,1,-1
380  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
381  C? Patrick Is this formula correct?  C? Patrick Is this formula correct?
382  cph Yes, but I rewrote it.  cph Yes, but I rewrote it.
383  cph Also, the KappaR? need the index and subscript k!  cph Also, the kappaR? need the index and subscript k!
384           kkey = (ikey-1)*Nr + k           kkey = (itdkey-1)*Nr + k
385  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
386    
387  C--       km1    Points to level above k (=k-1)  C--       km1    Points to level above k (=k-1)
# Line 521  C--       kDown  Cycles through 2,1 to p Line 397  C--       kDown  Cycles through 2,1 to p
397            jMin = 1-OLy            jMin = 1-OLy
398            jMax = sNy+OLy            jMax = sNy+OLy
399    
400              IF (k.EQ.Nr) THEN
401               DO j=1-Oly,sNy+Oly
402                DO i=1-Olx,sNx+Olx
403                 rTransKp1(i,j) = 0. _d 0
404                ENDDO
405               ENDDO
406              ELSE
407               DO j=1-Oly,sNy+Oly
408                DO i=1-Olx,sNx+Olx
409                 rTransKp1(i,j) = rTrans(i,j)
410                ENDDO
411               ENDDO
412              ENDIF
413    #ifdef ALLOW_AUTODIFF_TAMC
414    CADJ STORE rTransKp1(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
415    #endif
416    
417  C--       Get temporary terms used by tendency routines  C--       Get temporary terms used by tendency routines
418            CALL CALC_COMMON_FACTORS (            CALL CALC_COMMON_FACTORS (
419       I         bi,bj,iMin,iMax,jMin,jMax,k,       I         bi,bj,iMin,iMax,jMin,jMax,k,
420       O         xA,yA,uTrans,vTrans,rTrans,maskUp,       O         xA,yA,uTrans,vTrans,rTrans,maskUp,
421       I         myThid)       I         myThid)
422    
423              IF (k.EQ.1) THEN
424    C- Surface interface :
425               DO j=1-Oly,sNy+Oly
426                DO i=1-Olx,sNx+Olx
427                 rTrans(i,j) = 0.
428                ENDDO
429               ENDDO
430              ELSE
431    C- Interior interface :
432               DO j=1-Oly,sNy+Oly
433                DO i=1-Olx,sNx+Olx
434                 rTrans(i,j) = rTrans(i,j)*maskC(i,j,k-1,bi,bj)
435                ENDDO
436               ENDDO
437              ENDIF
438    
439  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
440    
441  C--   Residual transp = Bolus transp + Eulerian transp  C--   Residual transp = Bolus transp + Eulerian transp
442            IF (useGMRedi) THEN            IF (useGMRedi) THEN
443              CALL GMREDI_CALC_UVFLOW(              CALL GMREDI_CALC_UVFLOW(
# Line 535  C--   Residual transp = Bolus transp + E Line 445  C--   Residual transp = Bolus transp + E
445              IF (K.GE.2) CALL GMREDI_CALC_WFLOW(              IF (K.GE.2) CALL GMREDI_CALC_WFLOW(
446       &                    rTrans, bi, bj, k, myThid)       &                    rTrans, bi, bj, k, myThid)
447            ENDIF            ENDIF
 #endif /* ALLOW_GMREDI */  
448    
449  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
450  CADJ STORE KappaRT(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE rTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
451  CADJ STORE KappaRS(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte  #ifdef GM_BOLUS_ADVEC
452    CADJ STORE uTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
453    CADJ STORE vTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
454    #endif
455  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
456    
457    #endif /* ALLOW_GMREDI */
458    
459  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL
460  C--      Calculate the total vertical diffusivity  C--      Calculate the total vertical diffusivity
461           CALL CALC_DIFFUSIVITY(            IF ( .NOT.implicitDiffusion ) THEN
462       I        bi,bj,iMin,iMax,jMin,jMax,k,              CALL CALC_DIFFUSIVITY(
463       I        maskUp,       I          bi,bj,iMin,iMax,jMin,jMax,k,
464       O        KappaRT,KappaRS,       I          maskUp,
465       I        myThid)       O          kappaRT,kappaRS,
466         I          myThid)
467              ENDIF
468    # ifdef ALLOW_AUTODIFF_TAMC
469    CADJ STORE kappaRT(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
470    CADJ STORE kappaRS(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
471    # endif /* ALLOW_AUTODIFF_TAMC */
472  #endif  #endif
473    
474            iMin = 1-OLx+2            iMin = 1-OLx+2
# Line 557  C--      Calculate the total vertical di Line 477  C--      Calculate the total vertical di
477            jMax = sNy+OLy-1            jMax = sNy+OLy-1
478    
479  C--      Calculate active tracer tendencies (gT,gS,...)  C--      Calculate active tracer tendencies (gT,gS,...)
480  C        and step forward storing result in gTnm1, gSnm1, etc.  C        and step forward storing result in gT, gS, etc.
481    C--
482    # ifdef ALLOW_AUTODIFF_TAMC
483    #  if ((defined NONLIN_FRSURF) || (defined ALLOW_DEPTH_CONTROL)) && (defined ALLOW_GMREDI)
484    CADJ STORE kux(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
485    CADJ STORE kvy(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
486    #   ifdef GM_EXTRA_DIAGONAL
487    CADJ STORE kuz(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
488    CADJ STORE kvz(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
489    #   endif
490    #  endif
491    # endif /* ALLOW_AUTODIFF_TAMC */
492           IF ( tempStepping ) THEN           IF ( tempStepping ) THEN
493    #ifdef ALLOW_AUTODIFF_TAMC
494    CADJ STORE gTnm1(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
495    # if (defined NONLIN_FRSURF) || (defined ALLOW_DEPTH_CONTROL)
496    CADJ STORE gT(:,:,k,bi,bj)    = comlev1_bibj_k, key=kkey, byte=isbyte
497    # endif
498    #endif /* ALLOW_AUTODIFF_TAMC */
499             CALL CALC_GT(             CALL CALC_GT(
500       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
501       I         xA,yA,uTrans,vTrans,rTrans,maskUp,       I         xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,
502       I         KappaRT,       I         kappaRT,
503       U         fVerT,       U         fVerT,
504       I         myTime,myIter,myThid)       I         myTime,myIter,myThid)
505    #ifdef ALLOW_ADAMSBASHFORTH_3
506              IF ( AdamsBashforth_T ) THEN
507               CALL TIMESTEP_TRACER(
508         I         bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,
509         I         gtNm(1-Olx,1-Oly,1,1,1,m2),
510         U         gT,
511         I         myIter, myThid)
512              ELSE
513    #endif
514             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
515       I         bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,       I         bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,
516       I         theta, gT,       I         theta,
517         U         gT,
518       I         myIter, myThid)       I         myIter, myThid)
519    #ifdef ALLOW_ADAMSBASHFORTH_3
520              ENDIF
521    #endif
522           ENDIF           ENDIF
523    
524           IF ( saltStepping ) THEN           IF ( saltStepping ) THEN
525    #ifdef ALLOW_AUTODIFF_TAMC
526    CADJ STORE gSnm1(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
527    # if (defined NONLIN_FRSURF) || (defined ALLOW_DEPTH_CONTROL)
528    CADJ STORE gS(:,:,k,bi,bj)    = comlev1_bibj_k, key=kkey, byte=isbyte
529    # endif
530    #endif /* ALLOW_AUTODIFF_TAMC */
531             CALL CALC_GS(             CALL CALC_GS(
532       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
533       I         xA,yA,uTrans,vTrans,rTrans,maskUp,       I         xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,
534       I         KappaRS,       I         kappaRS,
535       U         fVerS,       U         fVerS,
536       I         myTime,myIter,myThid)       I         myTime,myIter,myThid)
537    #ifdef ALLOW_ADAMSBASHFORTH_3
538              IF ( AdamsBashforth_S ) THEN
539             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
540       I         bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,       I         bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,
541       I         salt, gS,       I         gsNm(1-Olx,1-Oly,1,1,1,m2),
542         U         gS,
543       I         myIter, myThid)       I         myIter, myThid)
544           ENDIF            ELSE
545  #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)  
546             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
547       I         bi,bj,iMin,iMax,jMin,jMax,k,tracerAdvScheme,       I         bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,
548       I         Tr1, gTr1,       I         salt,
549       I         myIter,myThid)       U         gS,
550           ENDIF       I         myIter, myThid)
551    #ifdef ALLOW_ADAMSBASHFORTH_3
552              ENDIF
553  #endif  #endif
554             ENDIF
555    
556  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
557           IF ( usePTRACERS ) THEN           IF ( usePTRACERS ) THEN
558             CALL PTRACERS_INTEGERATE(             IF ( .NOT.implicitDiffusion ) THEN
559                 CALL PTRACERS_CALC_DIFF(
560         I            bi,bj,iMin,iMax,jMin,jMax,k,
561         I            maskUp,
562         O            kappaRTr,
563         I            myThid)
564               ENDIF
565    # ifdef ALLOW_AUTODIFF_TAMC
566    CADJ STORE kappaRTr(:,:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
567    # endif /* ALLOW_AUTODIFF_TAMC */
568               CALL PTRACERS_INTEGRATE(
569       I         bi,bj,k,       I         bi,bj,k,
570       I         xA,yA,uTrans,vTrans,rTrans,maskUp,       I         xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,
571       X         KappaRS,       U         fVerP,
572         I         kappaRTr,
573       I         myIter,myTime,myThid)       I         myIter,myTime,myThid)
574           ENDIF           ENDIF
575  #endif /* ALLOW_PTRACERS */  #endif /* ALLOW_PTRACERS */
# Line 614  C--      Apply open boundary conditions Line 582  C--      Apply open boundary conditions
582  #endif   /* ALLOW_OBCS */  #endif   /* ALLOW_OBCS */
583    
584  C--      Freeze water  C--      Freeze water
585           IF (allowFreezing) THEN  C  this bit of code is left here for backward compatibility.
586    C  freezing at surface level has been moved to FORWARD_STEP
587             IF ( useOldFreezing .AND. .NOT. useSEAICE
588         &       .AND. .NOT.(useThSIce.AND.k.EQ.1) ) THEN
589  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
590  CADJ STORE gT(:,:,k,bi,bj) = comlev1_bibj_k  CADJ STORE gT(:,:,k,bi,bj) = comlev1_bibj_k
591  CADJ &   , key = kkey, byte = isbyte  CADJ &   , key = kkey, byte = isbyte
592  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
593              CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )              CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )
594           END IF           ENDIF
595    
596  C--     end of thermodynamic k loop (Nr:1)  C--     end of thermodynamic k loop (Nr:1)
597          ENDDO          ENDDO
598    
599    C       All explicit advection/diffusion/sources should now be
600  #ifdef ALLOW_AUTODIFF_TAMC  C       done. The updated tracer field is in gPtr. Accumalate
601  C? Patrick? What about this one?  C       explicit tendency and also reset gPtr to initial tracer
602  cph Keys iikey and idkey dont seem to be needed  C       field for implicit matrix calculation
603  cph since storing occurs on different tape for each  
604  cph impldiff call anyways.  #ifdef ALLOW_MATRIX
605  cph Thus, common block comlev1_impl isnt needed either.          IF (useMATRIX)
606  cph Storing below needed in the case useGMREDI.       &    CALL MATRIX_STORE_TENDENCY_EXP(bi,bj, myTime,myIter,myThid)
607          iikey = (ikey-1)*maximpl  #endif
608  #endif /* ALLOW_AUTODIFF_TAMC */  
609            iMin = 1
610  C--     Implicit diffusion          iMax = sNx
611          IF (implicitDiffusion) THEN          jMin = 1
612            jMax = sNy
613           IF (tempStepping) THEN  
614    C--     Implicit vertical advection & diffusion
615            IF ( tempStepping .AND. implicitDiffusion ) THEN
616              CALL CALC_3D_DIFFUSIVITY(
617         I         bi,bj,iMin,iMax,jMin,jMax,
618         I         GAD_TEMPERATURE, useGMredi, useKPP,
619         O         kappaRk,
620         I         myThid)
621            ENDIF
622    #ifdef INCLUDE_IMPLVERTADV_CODE
623            IF ( tempImplVertAdv ) THEN
624              CALL GAD_IMPLICIT_R(
625         I         tempImplVertAdv, tempAdvScheme, GAD_TEMPERATURE,
626         I         kappaRk, wVel, theta,
627         U         gT,
628         I         bi, bj, myTime, myIter, myThid )
629            ELSEIF ( tempStepping .AND. implicitDiffusion ) THEN
630    #else /* INCLUDE_IMPLVERTADV_CODE */
631            IF     ( tempStepping .AND. implicitDiffusion ) THEN
632    #endif /* INCLUDE_IMPLVERTADV_CODE */
633  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
634              idkey = iikey + 1  CADJ STORE kappaRk(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
635  CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
636  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
637              CALL IMPLDIFF(            CALL IMPLDIFF(
638       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
639       I         deltaTtracer, KappaRT, recip_HFacC,       I         GAD_TEMPERATURE, kappaRk, recip_hFacC,
640       U         gT,       U         gT,
641       I         myThid )       I         myThid )
642            ENDIF
643    
644    #ifdef ALLOW_TIMEAVE
645            useVariableK = useKPP .OR. usePP81 .OR. useMY82 .OR. useGGL90
646         &       .OR. useGMredi .OR. ivdc_kappa.NE.0.
647            IF (taveFreq.GT.0. .AND. useVariableK ) THEN
648             IF (implicitDiffusion) THEN
649               CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, gT, kappaRk,
650         I                        Nr, 3, deltaTclock, bi, bj, myThid)
651    c        ELSE
652    c          CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, theta, kappaRT,
653    c    I                        Nr, 3, deltaTclock, bi, bj, myThid)
654           ENDIF           ENDIF
655            ENDIF
656    #endif /* ALLOW_TIMEAVE */
657    
658            IF ( saltStepping .AND. implicitDiffusion ) THEN
659              CALL CALC_3D_DIFFUSIVITY(
660         I         bi,bj,iMin,iMax,jMin,jMax,
661         I         GAD_SALINITY, useGMredi, useKPP,
662         O         kappaRk,
663         I         myThid)
664            ENDIF
665    
666           IF (saltStepping) THEN  #ifdef INCLUDE_IMPLVERTADV_CODE
667            IF ( saltImplVertAdv ) THEN
668              CALL GAD_IMPLICIT_R(
669         I         saltImplVertAdv, saltAdvScheme, GAD_SALINITY,
670         I         kappaRk, wVel, salt,
671         U         gS,
672         I         bi, bj, myTime, myIter, myThid )
673            ELSEIF ( saltStepping .AND. implicitDiffusion ) THEN
674    #else /* INCLUDE_IMPLVERTADV_CODE */
675            IF     ( saltStepping .AND. implicitDiffusion ) THEN
676    #endif /* INCLUDE_IMPLVERTADV_CODE */
677  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
678           idkey = iikey + 2  CADJ STORE kappaRk(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
679  CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
680  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
681              CALL IMPLDIFF(            CALL IMPLDIFF(
682       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
683       I         deltaTtracer, KappaRS, recip_HFacC,       I         GAD_SALINITY, kappaRk, recip_hFacC,
684       U         gS,       U         gS,
685       I         myThid )       I         myThid )
686           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  
687    
688  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
689  C Vertical diffusion (implicit) for passive tracers          IF     ( usePTRACERS ) THEN
690           IF ( usePTRACERS ) THEN  C--     Vertical advection/diffusion (implicit) for passive tracers
691             CALL PTRACERS_IMPLDIFF( bi,bj,KappaRS,myThid )             CALL PTRACERS_IMPLICIT(
692           ENDIF       U                             kappaRk,
693         I                             bi, bj, myTime, myIter, myThid )
694            ENDIF
695  #endif /* ALLOW_PTRACERS */  #endif /* ALLOW_PTRACERS */
696    
697  #ifdef   ALLOW_OBCS  #ifdef   ALLOW_OBCS
698  C--      Apply open boundary conditions  C--      Apply open boundary conditions
699           IF (useOBCS) THEN          IF ( ( implicitDiffusion
700         &    .OR. tempImplVertAdv
701         &    .OR. saltImplVertAdv
702         &       ) .AND. useOBCS     ) THEN
703             DO K=1,Nr             DO K=1,Nr
704               CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )               CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )
705             ENDDO             ENDDO
          END IF  
 #endif   /* ALLOW_OBCS */  
   
 C--     End If implicitDiffusion  
706          ENDIF          ENDIF
707    #endif   /* ALLOW_OBCS */
708    
709  #endif /* SINGLE_LAYER_MODE */  #endif /* SINGLE_LAYER_MODE */
710    
711  Ccs-  C--   end bi,bj loops.
712         ENDDO         ENDDO
713        ENDDO        ENDDO
714    
715  #ifdef ALLOW_AIM  #ifdef ALLOW_DEBUG
       IF ( useAIM ) THEN  
        CALL AIM_AIM2DYN_EXCHANGES( myTime, myThid )  
       ENDIF  
 #endif /* ALLOW_AIM */  
       IF ( staggerTimeStep ) THEN  
        IF ( useAIM .OR. useCubedSphereExchange ) THEN  
          IF (tempStepping) _EXCH_XYZ_R8(gT,myThid)  
          IF (saltStepping) _EXCH_XYZ_R8(gS,myThid)  
        ELSEIF ( useGMRedi .AND. Oly.LT.4 ) THEN  
 c        .AND. GM_AdvForm .AND. .NOT.GM_AdvSeparate ) THEN  
          IF (tempMultiDimAdvec) _EXCH_XYZ_R8(gT,myThid)  
          IF (saltMultiDimAdvec) _EXCH_XYZ_R8(gS,myThid)  
        ENDIF    
       ENDIF    
   
 #ifndef DISABLE_DEBUGMODE  
716        If (debugMode) THEN        If (debugMode) THEN
717         CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (THERMODYNAMICS)',myThid)
718         CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (THERMODYNAMICS)',myThid)
719         CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (THERMODYNAMICS)',myThid)
720         CALL DEBUG_STATS_RL(Nr,theta,'Theta (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,theta,'Theta (THERMODYNAMICS)',myThid)
721         CALL DEBUG_STATS_RL(Nr,salt,'Salt (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,salt,'Salt (THERMODYNAMICS)',myThid)
722         CALL DEBUG_STATS_RL(Nr,Gt,'Gt (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,gT,'Gt (THERMODYNAMICS)',myThid)
723         CALL DEBUG_STATS_RL(Nr,Gs,'Gs (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,gS,'Gs (THERMODYNAMICS)',myThid)
724         CALL DEBUG_STATS_RL(Nr,GtNm1,'GtNm1 (THERMODYNAMICS)',myThid)  #ifndef ALLOW_ADAMSBASHFORTH_3
725         CALL DEBUG_STATS_RL(Nr,GsNm1,'GsNm1 (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,gtNm1,'GtNm1 (THERMODYNAMICS)',myThid)
726           CALL DEBUG_STATS_RL(Nr,gsNm1,'GsNm1 (THERMODYNAMICS)',myThid)
727    #endif
728  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
729         IF ( usePTRACERS ) THEN         IF ( usePTRACERS ) THEN
730           CALL PTRACERS_DEBUG(myThid)           CALL PTRACERS_DEBUG(myThid)
731         ENDIF         ENDIF
732  #endif /* ALLOW_PTRACERS */  #endif /* ALLOW_PTRACERS */
733        ENDIF        ENDIF
734    #endif /* ALLOW_DEBUG */
735    
736    #ifdef ALLOW_DEBUG
737             IF ( debugLevel .GE. debLevB )
738         &    CALL DEBUG_LEAVE('THERMODYNAMICS',myThid)
739  #endif  #endif
740    
741    #endif /* ALLOW_GENERIC_ADVDIFF */
742    
743        RETURN        RETURN
744        END        END
Legend:
Removed from v.1.24  
changed lines
  Added in v.1.104
  ViewVC Help
Powered by ViewVC 1.1.22