/[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.10 by adcroft, Thu Sep 27 18:06:43 2001 UTC revision 1.50 by jmc, Tue Oct 7 04:31:30 2003 UTC
# Line 2  C $Header$ Line 2  C $Header$
2  C $Name$  C $Name$
3    
4  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
5    #ifdef ALLOW_AUTODIFF_TAMC
6    # ifdef ALLOW_GMREDI
7    #  include "GMREDI_OPTIONS.h"
8    # endif
9    # ifdef ALLOW_KPP
10    #  include "KPP_OPTIONS.h"
11    # endif
12    #ifdef ALLOW_PTRACERS
13    # include "PTRACERS_OPTIONS.h"
14    #endif
15    #endif /* ALLOW_AUTODIFF_TAMC */
16    
17  CBOP  CBOP
18  C     !ROUTINE: THERMODYNAMICS  C     !ROUTINE: THERMODYNAMICS
# Line 71  C     == Global variables === Line 82  C     == Global variables ===
82  #ifdef ALLOW_PASSIVE_TRACER  #ifdef ALLOW_PASSIVE_TRACER
83  #include "TR1.h"  #include "TR1.h"
84  #endif  #endif
85    #ifdef ALLOW_PTRACERS
86    #include "PTRACERS.h"
87    #endif
88    #ifdef ALLOW_TIMEAVE
89    #include "TIMEAVE_STATV.h"
90    #endif
91    
92  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
93  # include "tamc.h"  # include "tamc.h"
94  # include "tamc_keys.h"  # include "tamc_keys.h"
95  # include "FFIELDS.h"  # include "FFIELDS.h"
96    # include "EOS.h"
97  # ifdef ALLOW_KPP  # ifdef ALLOW_KPP
98  #  include "KPP.h"  #  include "KPP.h"
99  # endif  # endif
# Line 82  C     == Global variables === Line 101  C     == Global variables ===
101  #  include "GMREDI.h"  #  include "GMREDI.h"
102  # endif  # endif
103  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
 #ifdef ALLOW_TIMEAVE  
 #include "TIMEAVE_STATV.h"  
 #endif  
104    
105  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
106  C     == Routine arguments ==  C     == Routine arguments ==
# Line 109  C                                      i Line 125  C                                      i
125  C                                      so we need an fVer for each  C                                      so we need an fVer for each
126  C                                      variable.  C                                      variable.
127  C     rhoK, rhoKM1   - Density at current level, and level above  C     rhoK, rhoKM1   - Density at current level, and level above
 C     phiHyd         - Hydrostatic part of the potential phiHydi.  
 C                      In z coords phiHydiHyd is the hydrostatic  
 C                      Potential (=pressure/rho0) anomaly  
 C                      In p coords phiHydiHyd is the geopotential  
 C                      surface height anomaly.  
128  C     phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean)  C     phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean)
129  C     phiSurfY             or geopotentiel (atmos) in X and Y direction  C     phiSurfY             or geopotentiel (atmos) in X and Y direction
130  C     KappaRT,       - Total diffusion in vertical for T and S.  C     KappaRT,       - Total diffusion in vertical for T and S.
131  C     KappaRS          (background + spatially varying, isopycnal term).  C     KappaRS          (background + spatially varying, isopycnal term).
132    C     useVariableK   = T when vertical diffusion is not constant
133  C     iMin, iMax     - Ranges and sub-block indices on which calculations  C     iMin, iMax     - Ranges and sub-block indices on which calculations
134  C     jMin, jMax       are applied.  C     jMin, jMax       are applied.
135  C     bi, bj  C     bi, bj
# Line 133  C                      index into fVerTe Line 145  C                      index into fVerTe
145        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
146        _RL fVerS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
147        _RL fVerTr1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerTr1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
       _RL phiHyd  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
148        _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
149        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
150        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 145  C                      index into fVerTe Line 156  C                      index into fVerTe
156        _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
157  C     This is currently used by IVDC and Diagnostics  C     This is currently used by IVDC and Diagnostics
158        _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
159          LOGICAL useVariableK
160        INTEGER iMin, iMax        INTEGER iMin, iMax
161        INTEGER jMin, jMax        INTEGER jMin, jMax
162        INTEGER bi, bj        INTEGER bi, bj
163        INTEGER i, j        INTEGER i, j
164        INTEGER k, km1, kup, kDown        INTEGER k, km1, kup, kDown
165          INTEGER iTracer
166    
 Cjmc : add for phiHyd output <- but not working if multi tile per CPU  
 c     CHARACTER*(MAX_LEN_MBUF) suff  
 c     LOGICAL  DIFFERENT_MULTIPLE  
 c     EXTERNAL DIFFERENT_MULTIPLE  
 Cjmc(end)  
167  CEOP  CEOP
168    
169    #ifndef DISABLE_DEBUGMODE
170             IF ( debugLevel .GE. debLevB )
171         &    CALL DEBUG_ENTER('FORWARD_STEP',myThid)
172    #endif
173    
174  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
175  C--   dummy statement to end declaration part  C--   dummy statement to end declaration part
176        ikey = 1        ikey = 1
177          itdkey = 1
178  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
179    
 C--   Set up work arrays with valid (i.e. not NaN) values  
 C     These inital values do not alter the numerical results. They  
 C     just ensure that all memory references are to valid floating  
 C     point numbers. This prevents spurious hardware signals due to  
 C     uninitialised but inert locations.  
       DO j=1-OLy,sNy+OLy  
        DO i=1-OLx,sNx+OLx  
         xA(i,j)      = 0. _d 0  
         yA(i,j)      = 0. _d 0  
         uTrans(i,j)  = 0. _d 0  
         vTrans(i,j)  = 0. _d 0  
         DO k=1,Nr  
          phiHyd(i,j,k)  = 0. _d 0  
          sigmaX(i,j,k) = 0. _d 0  
          sigmaY(i,j,k) = 0. _d 0  
          sigmaR(i,j,k) = 0. _d 0  
         ENDDO  
         rhoKM1 (i,j) = 0. _d 0  
         rhok   (i,j) = 0. _d 0  
         phiSurfX(i,j) = 0. _d 0  
         phiSurfY(i,j) = 0. _d 0  
        ENDDO  
       ENDDO  
   
   
180  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
181  C--   HPF directive to help TAMC  C--   HPF directive to help TAMC
182  CHPF$ INDEPENDENT  CHPF$ INDEPENDENT
# Line 198  CHPF$ INDEPENDENT Line 187  CHPF$ INDEPENDENT
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, NEW (rTrans,fVerT,fVerS  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS
190  CHPF$&                  ,phiHyd,utrans,vtrans,xA,yA  CHPF$&                  ,utrans,vtrans,xA,yA
191  CHPF$&                  ,KappaRT,KappaRS  CHPF$&                  ,KappaRT,KappaRS
192  CHPF$&                  )  CHPF$&                  )
193  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
# Line 208  CHPF$&                  ) Line 197  CHPF$&                  )
197  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
198            act1 = bi - myBxLo(myThid)            act1 = bi - myBxLo(myThid)
199            max1 = myBxHi(myThid) - myBxLo(myThid) + 1            max1 = myBxHi(myThid) - myBxLo(myThid) + 1
   
200            act2 = bj - myByLo(myThid)            act2 = bj - myByLo(myThid)
201            max2 = myByHi(myThid) - myByLo(myThid) + 1            max2 = myByHi(myThid) - myByLo(myThid) + 1
   
202            act3 = myThid - 1            act3 = myThid - 1
203            max3 = nTx*nTy            max3 = nTx*nTy
   
204            act4 = ikey_dynamics - 1            act4 = ikey_dynamics - 1
205              itdkey = (act1 + 1) + act2*max1
           ikey = (act1 + 1) + act2*max1  
206       &                      + act3*max1*max2       &                      + act3*max1*max2
207       &                      + act4*max1*max2*max3       &                      + act4*max1*max2*max3
208  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
209    
210  C--     Set up work arrays that need valid initial values  C--   Set up work arrays with valid (i.e. not NaN) values
211    C     These inital values do not alter the numerical results. They
212    C     just ensure that all memory references are to valid floating
213    C     point numbers. This prevents spurious hardware signals due to
214    C     uninitialised but inert locations.
215    
216          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
217           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
218              xA(i,j)        = 0. _d 0
219              yA(i,j)        = 0. _d 0
220              uTrans(i,j)    = 0. _d 0
221              vTrans(i,j)    = 0. _d 0
222              rhok   (i,j)   = 0. _d 0
223              rhoKM1 (i,j)   = 0. _d 0
224              phiSurfX(i,j)  = 0. _d 0
225              phiSurfY(i,j)  = 0. _d 0
226            rTrans (i,j)   = 0. _d 0            rTrans (i,j)   = 0. _d 0
227            fVerT  (i,j,1) = 0. _d 0            fVerT  (i,j,1) = 0. _d 0
228            fVerT  (i,j,2) = 0. _d 0            fVerT  (i,j,2) = 0. _d 0
# Line 239  C--     Set up work arrays that need val Line 237  C--     Set up work arrays that need val
237           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
238            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
239  C This is currently also used by IVDC and Diagnostics  C This is currently also used by IVDC and Diagnostics
240               sigmaX(i,j,k) = 0. _d 0
241               sigmaY(i,j,k) = 0. _d 0
242               sigmaR(i,j,k) = 0. _d 0
243             ConvectCount(i,j,k) = 0.             ConvectCount(i,j,k) = 0.
244             KappaRT(i,j,k) = 0. _d 0             KappaRT(i,j,k)    = 0. _d 0
245             KappaRS(i,j,k) = 0. _d 0             KappaRS(i,j,k)    = 0. _d 0
246  #ifdef ALLOW_AUTODIFF_TAMC  C- tracer tendency needs to be set to zero (moved here from gad_calc_rhs):
247             gT(i,j,k,bi,bj) = 0. _d 0             gT(i,j,k,bi,bj)   = 0. _d 0
248             gS(i,j,k,bi,bj) = 0. _d 0             gS(i,j,k,bi,bj)   = 0. _d 0
249  #ifdef ALLOW_PASSIVE_TRACER  # ifdef ALLOW_PASSIVE_TRACER
250             gTr1(i,j,k,bi,bj) = 0. _d 0             gTr1(i,j,k,bi,bj) = 0. _d 0
251  #endif  # endif
252  #endif  # ifdef ALLOW_PTRACERS
253               DO iTracer=1,PTRACERS_numInUse
254                gPTr(i,j,k,bi,bj,itracer) = 0. _d 0
255               ENDDO
256    # endif
257    #ifdef ALLOW_AUTODIFF_TAMC
258    cph all the following init. are necessary for TAF
259    cph although some of these are re-initialised later.
260    # ifdef ALLOW_GMREDI
261               Kwx(i,j,k,bi,bj)  = 0. _d 0
262               Kwy(i,j,k,bi,bj)  = 0. _d 0
263               Kwz(i,j,k,bi,bj)  = 0. _d 0
264    #  ifdef GM_NON_UNITY_DIAGONAL
265               Kux(i,j,k,bi,bj)  = 0. _d 0
266               Kvy(i,j,k,bi,bj)  = 0. _d 0
267    #  endif
268    #  ifdef GM_EXTRA_DIAGONAL
269               Kuz(i,j,k,bi,bj)  = 0. _d 0
270               Kvz(i,j,k,bi,bj)  = 0. _d 0
271    #  endif
272    #  ifdef GM_BOLUS_ADVEC
273               GM_PsiX(i,j,k,bi,bj)  = 0. _d 0
274               GM_PsiY(i,j,k,bi,bj)  = 0. _d 0
275    #  endif
276    #  ifdef GM_VISBECK_VARIABLE_K
277               VisbeckK(i,j,bi,bj)   = 0. _d 0
278    #  endif
279    # endif /* ALLOW_GMREDI */
280    #endif /* ALLOW_AUTODIFF_TAMC */
281            ENDDO            ENDDO
282           ENDDO           ENDDO
283          ENDDO          ENDDO
284    
285          iMin = 1-OLx+1          iMin = 1-OLx
286          iMax = sNx+OLx          iMax = sNx+OLx
287          jMin = 1-OLy+1          jMin = 1-OLy
288          jMax = sNy+OLy          jMax = sNy+OLy
289    
   
290  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
291  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
292  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
293    CADJ STORE totphihyd
294    CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
295    #ifdef ALLOW_KPP
296    CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
297    CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
298    #endif
299  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
300    
301    #ifndef DISABLE_DEBUGMODE
302            IF ( debugLevel .GE. debLevB )
303         &    CALL DEBUG_MSG('ENTERING UPWARD K LOOP',myThid)
304    #endif
305    
306  C--     Start of diagnostic loop  C--     Start of diagnostic loop
307          DO k=Nr,1,-1          DO k=Nr,1,-1
308    
# Line 272  C? Patrick, is this formula correct now Line 311  C? Patrick, is this formula correct now
311  C? Do we still need this?  C? Do we still need this?
312  cph kkey formula corrected.  cph kkey formula corrected.
313  cph Needed for rhok, rhokm1, in the case useGMREDI.  cph Needed for rhok, rhokm1, in the case useGMREDI.
314           kkey = (ikey-1)*Nr + k           kkey = (itdkey-1)*Nr + k
 CADJ STORE rhokm1(:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte  
 CADJ STORE rhok  (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte  
315  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
316    
317  C--       Integrate continuity vertically for vertical velocity  C--       Integrate continuity vertically for vertical velocity
318            CALL INTEGRATE_FOR_W(  c         CALL INTEGRATE_FOR_W(
319       I                         bi, bj, k, uVel, vVel,  c    I                         bi, bj, k, uVel, vVel,
320       O                         wVel,  c    O                         wVel,
321       I                         myThid )  c    I                         myThid )
322    
323  #ifdef    ALLOW_OBCS  #ifdef    ALLOW_OBCS
324  #ifdef    ALLOW_NONHYDROSTATIC  #ifdef    ALLOW_NONHYDROSTATIC
325  C--       Apply OBC to W if in N-H mode  C--       Apply OBC to W if in N-H mode
326            IF (useOBCS.AND.nonHydrostatic) THEN  c         IF (useOBCS.AND.nonHydrostatic) THEN
327              CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )  c           CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )
328            ENDIF  c         ENDIF
329  #endif    /* ALLOW_NONHYDROSTATIC */  #endif    /* ALLOW_NONHYDROSTATIC */
330  #endif    /* ALLOW_OBCS */  #endif    /* ALLOW_OBCS */
331    
332    C--       Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h
333    C--       MOST of THERMODYNAMICS will be disabled
334    #ifndef SINGLE_LAYER_MODE
335    
336  C--       Calculate gradients of potential density for isoneutral  C--       Calculate gradients of potential density for isoneutral
337  C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)  C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)
338  c         IF ( k.GT.1 .AND. (useGMRedi.OR.ivdc_kappa.NE.0.) ) THEN  c         IF ( k.GT.1 .AND. (useGMRedi.OR.ivdc_kappa.NE.0.) ) THEN
339            IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.) ) THEN            IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.) ) THEN
340    #ifndef DISABLE_DEBUGMODE
341                IF ( debugLevel .GE. debLevB )
342         &       CALL DEBUG_CALL('FIND_RHO',myThid)
343    #endif
344  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
345  CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
346  CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
347  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
348              CALL FIND_RHO(              CALL FIND_RHO(
349       I        bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,       I        bi, bj, iMin, iMax, jMin, jMax, k, k,
350       I        theta, salt,       I        theta, salt,
351       O        rhoK,       O        rhoK,
352       I        myThid )       I        myThid )
353    
354              IF (k.GT.1) THEN              IF (k.GT.1) THEN
355  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
356  CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
357  CADJ STORE salt (:,:,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
358  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
359               CALL FIND_RHO(               CALL FIND_RHO(
360       I        bi, bj, iMin, iMax, jMin, jMax, k-1, k, eosType,       I        bi, bj, iMin, iMax, jMin, jMax, k-1, k,
361       I        theta, salt,       I        theta, salt,
362       O        rhoKm1,       O        rhoKm1,
363       I        myThid )       I        myThid )
364              ENDIF              ENDIF
365    #ifndef DISABLE_DEBUGMODE
366                IF ( debugLevel .GE. debLevB )
367         &       CALL DEBUG_CALL('GRAD_SIGMA',myThid)
368    #endif
369              CALL GRAD_SIGMA(              CALL GRAD_SIGMA(
370       I             bi, bj, iMin, iMax, jMin, jMax, k,       I             bi, bj, iMin, iMax, jMin, jMax, k,
371       I             rhoK, rhoKm1, rhoK,       I             rhoK, rhoKm1, rhoK,
# Line 323  CADJ STORE salt (:,:,k-1,bi,bj) = comlev Line 373  CADJ STORE salt (:,:,k-1,bi,bj) = comlev
373       I             myThid )       I             myThid )
374            ENDIF            ENDIF
375    
376    #ifdef ALLOW_AUTODIFF_TAMC
377    CADJ STORE rhok   (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte
378    CADJ STORE rhokm1 (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte
379    #endif /* ALLOW_AUTODIFF_TAMC */
380  C--       Implicit Vertical Diffusion for Convection  C--       Implicit Vertical Diffusion for Convection
381  c ==> should use sigmaR !!!  c ==> should use sigmaR !!!
382            IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN            IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN
383    #ifndef DISABLE_DEBUGMODE
384                IF ( debugLevel .GE. debLevB )
385         &       CALL DEBUG_CALL('CALC_IVDC',myThid)
386    #endif
387              CALL CALC_IVDC(              CALL CALC_IVDC(
388       I        bi, bj, iMin, iMax, jMin, jMax, k,       I        bi, bj, iMin, iMax, jMin, jMax, k,
389       I        rhoKm1, rhoK,       I        rhoKm1, rhoK,
# Line 333  c ==> should use sigmaR !!! Line 391  c ==> should use sigmaR !!!
391       I        myTime, myIter, myThid)       I        myTime, myIter, myThid)
392            ENDIF            ENDIF
393    
394    #endif /* SINGLE_LAYER_MODE */
395    
396  C--     end of diagnostic k loop (Nr:1)  C--     end of diagnostic k loop (Nr:1)
397          ENDDO          ENDDO
398    
399  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
400  cph avoids recomputation of integrate_for_w  cph avoids recomputation of integrate_for_w
401  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
402  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
403    
404  #ifdef  ALLOW_OBCS  #ifdef  ALLOW_OBCS
405  C--     Calculate future values on open boundaries  C--     Calculate future values on open boundaries
406          IF (useOBCS) THEN          IF (useOBCS) THEN
407            CALL OBCS_CALC( bi, bj, myTime+deltaT,  #ifndef DISABLE_DEBUGMODE
408              IF ( debugLevel .GE. debLevB )
409         &     CALL DEBUG_CALL('OBCS_CALC',myThid)
410    #endif
411              CALL OBCS_CALC( bi, bj, myTime+deltaT, myIter+1,
412       I            uVel, vVel, wVel, theta, salt,       I            uVel, vVel, wVel, theta, salt,
413       I            myThid )       I            myThid )
414          ENDIF          ENDIF
415  #endif  /* ALLOW_OBCS */  #endif  /* ALLOW_OBCS */
416    
417    
418    #ifdef ALLOW_THERM_SEAICE
419           IF (useThermSeaIce) THEN
420    #ifndef DISABLE_DEBUGMODE
421            IF ( debugLevel .GE. debLevB )
422         &    CALL DEBUG_CALL('ICE_FORCING',myThid)
423    #endif
424  C--     Determines forcing terms based on external fields  C--     Determines forcing terms based on external fields
425  C       relaxation terms, etc.  C       including effects from ice
426          CALL EXTERNAL_FORCING_SURF(          CALL ICE_FORCING(
427       I             bi, bj, iMin, iMax, jMin, jMax,       I             bi, bj, iMin, iMax, jMin, jMax,
428       I             myThid )       I             myThid )
429           ELSE
430    #endif /* ALLOW_THERM_SEAICE */
431    
432    C--     Determines forcing terms based on external fields
433    C       relaxation terms, etc.
434    #ifndef DISABLE_DEBUGMODE
435            IF ( debugLevel .GE. debLevB )
436         &    CALL DEBUG_CALL('EXTERNAL_FORCING_SURF',myThid)
437    #endif
438            CALL EXTERNAL_FORCING_SURF(
439         I             bi, bj, iMin, iMax, jMin, jMax,
440         I             myTime, myIter, myThid )
441    
442    #ifdef ALLOW_THERM_SEAICE
443    C--    end of if/else block useThermSeaIce --
444           ENDIF
445    #endif /* ALLOW_THERM_SEAICE */
446    
447  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
448  cph needed for KPP  cph needed for KPP
449  CADJ STORE surfacetendencyU(:,:,bi,bj)  CADJ STORE surfacetendencyU(:,:,bi,bj)
450  CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
451  CADJ STORE surfacetendencyV(:,:,bi,bj)  CADJ STORE surfacetendencyV(:,:,bi,bj)
452  CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
453  CADJ STORE surfacetendencyS(:,:,bi,bj)  CADJ STORE surfacetendencyS(:,:,bi,bj)
454  CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
455  CADJ STORE surfacetendencyT(:,:,bi,bj)  CADJ STORE surfacetendencyT(:,:,bi,bj)
456  CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
457  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
458    
459    C--     Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h
460    C--     MOST of THERMODYNAMICS will be disabled
461    #ifndef SINGLE_LAYER_MODE
462    
463  #ifdef  ALLOW_GMREDI  #ifdef  ALLOW_GMREDI
464    
465  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
466  CADJ STORE sigmaX(:,:,:) = comlev1, key=ikey, byte=isbyte  cph storing here is needed only for one GMREDI_OPTIONS:
467  CADJ STORE sigmaY(:,:,:) = comlev1, key=ikey, byte=isbyte  cph define GM_BOLUS_ADVEC
468  CADJ STORE sigmaR(:,:,:) = comlev1, key=ikey, byte=isbyte  cph but I've avoided the #ifdef for now, in case more things change
469    CADJ STORE sigmaX(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte
470    CADJ STORE sigmaY(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte
471    CADJ STORE sigmaR(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte
472  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
473    
474  C--     Calculate iso-neutral slopes for the GM/Redi parameterisation  C--     Calculate iso-neutral slopes for the GM/Redi parameterisation
475          IF (useGMRedi) THEN          IF (useGMRedi) THEN
476            DO k=1,Nr  #ifndef DISABLE_DEBUGMODE
477              CALL GMREDI_CALC_TENSOR(            IF ( debugLevel .GE. debLevB )
478       I             bi, bj, iMin, iMax, jMin, jMax, k,       &     CALL DEBUG_CALL('GMREDI_CALC_TENSOR',myThid)
479    #endif
480              CALL GMREDI_CALC_TENSOR(
481         I             bi, bj, iMin, iMax, jMin, jMax,
482       I             sigmaX, sigmaY, sigmaR,       I             sigmaX, sigmaY, sigmaR,
483       I             myThid )       I             myThid )
           ENDDO  
484  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
485          ELSE          ELSE
486            DO k=1, Nr            CALL GMREDI_CALC_TENSOR_DUMMY(
487              CALL GMREDI_CALC_TENSOR_DUMMY(       I             bi, bj, iMin, iMax, jMin, jMax,
      I             bi, bj, iMin, iMax, jMin, jMax, k,  
488       I             sigmaX, sigmaY, sigmaR,       I             sigmaX, sigmaY, sigmaR,
489       I             myThid )       I             myThid )
           ENDDO  
490  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
491          ENDIF          ENDIF
492    
493  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
494  CADJ STORE Kwx(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE Kwx(:,:,:,bi,bj)   = comlev1_bibj, key=itdkey, byte=isbyte
495  CADJ STORE Kwy(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE Kwy(:,:,:,bi,bj)   = comlev1_bibj, key=itdkey, byte=isbyte
496  CADJ STORE Kwz(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE Kwz(:,:,:,bi,bj)   = comlev1_bibj, key=itdkey, byte=isbyte
497  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
498    
499  #endif  /* ALLOW_GMREDI */  #endif  /* ALLOW_GMREDI */
# Line 404  CADJ STORE Kwz(:,:,:,bi,bj)   = comlev1_ Line 501  CADJ STORE Kwz(:,:,:,bi,bj)   = comlev1_
501  #ifdef  ALLOW_KPP  #ifdef  ALLOW_KPP
502  C--     Compute KPP mixing coefficients  C--     Compute KPP mixing coefficients
503          IF (useKPP) THEN          IF (useKPP) THEN
504    #ifndef DISABLE_DEBUGMODE
505              IF ( debugLevel .GE. debLevB )
506         &     CALL DEBUG_CALL('KPP_CALC',myThid)
507    #endif
508            CALL KPP_CALC(            CALL KPP_CALC(
509       I                  bi, bj, myTime, myThid )       I                  bi, bj, myTime, myThid )
510  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 415  C--     Compute KPP mixing coefficients Line 516  C--     Compute KPP mixing coefficients
516    
517  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
518  CADJ STORE KPPghat   (:,:,:,bi,bj)  CADJ STORE KPPghat   (:,:,:,bi,bj)
 CADJ &   , KPPviscAz (:,:,:,bi,bj)  
519  CADJ &   , KPPdiffKzT(:,:,:,bi,bj)  CADJ &   , KPPdiffKzT(:,:,:,bi,bj)
520  CADJ &   , KPPdiffKzS(:,:,:,bi,bj)  CADJ &   , KPPdiffKzS(:,:,:,bi,bj)
521  CADJ &   , KPPfrac   (:,:  ,bi,bj)  CADJ &   , KPPfrac   (:,:  ,bi,bj)
522  CADJ &                 = comlev1_bibj, key=ikey, byte=isbyte  CADJ &                 = comlev1_bibj, key=itdkey, byte=isbyte
523  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
524    
525  #endif  /* ALLOW_KPP */  #endif  /* ALLOW_KPP */
526    
527  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
528  CADJ STORE KappaRT(:,:,:)     = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE KappaRT(:,:,:)     = comlev1_bibj, key=itdkey, byte=isbyte
529  CADJ STORE KappaRS(:,:,:)     = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE KappaRS(:,:,:)     = comlev1_bibj, key=itdkey, byte=isbyte
530  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
531  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
532  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
533  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
534  #ifdef ALLOW_PASSIVE_TRACER  #ifdef ALLOW_PASSIVE_TRACER
535  CADJ STORE tr1  (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE tr1  (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
536    #endif
537    #ifdef ALLOW_PTRACERS
538    cph-- moved to forward_step to avoid key computation
539    cphCADJ STORE ptracer(:,:,:,bi,bj,itracer) = comlev1_bibj,
540    cphCADJ &                              key=itdkey, byte=isbyte
541  #endif  #endif
542  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
543    
544  #ifdef ALLOW_AIM  #ifdef ALLOW_AIM
545  C       AIM - atmospheric intermediate model, physics package code.  C       AIM - atmospheric intermediate model, physics package code.
 C note(jmc) : phiHyd=0 at this point but is not really used in Molteni Physics  
546          IF ( useAIM ) THEN          IF ( useAIM ) THEN
547           CALL TIMER_START('AIM_DO_ATMOS_PHYS      [DYNAMICS]', myThid)  #ifndef DISABLE_DEBUGMODE
548           CALL AIM_DO_ATMOS_PHYSICS( phiHyd, bi, bj, myTime, myThid )            IF ( debugLevel .GE. debLevB )
549           CALL TIMER_STOP ('AIM_DO_ATMOS_PHYS      [DYNAMICS]', myThid)       &     CALL DEBUG_CALL('AIM_DO_PHYSICS',myThid)
550    #endif
551             CALL TIMER_START('AIM_DO_PHYSICS   [THERMODYNAMICS]', myThid)
552             CALL AIM_DO_PHYSICS( bi, bj, myTime, myIter, myThid )
553             CALL TIMER_STOP( 'AIM_DO_PHYSICS   [THERMODYNAMICS]', myThid)
554          ENDIF          ENDIF
555  #endif /* ALLOW_AIM */  #endif /* ALLOW_AIM */
556    
557    #ifndef DISABLE_MULTIDIM_ADVECTION
558  C--     Some advection schemes are better calculated using a multi-dimensional  C--     Some advection schemes are better calculated using a multi-dimensional
559  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.
560          IF (multiDimAdvection) THEN  C
561           IF (tempStepping .AND.  C The CPP flag DISABLE_MULTIDIM_ADVECTION is currently unset in GAD_OPTIONS.h
562       &       tempAdvScheme.NE.ENUM_CENTERED_2ND .AND.  C The default is to use multi-dimensinal advection for non-linear advection
563       &       tempAdvScheme.NE.ENUM_UPWIND_3RD .AND.  C schemes. However, for the sake of efficiency of the adjoint it is necessary
564       &       tempAdvScheme.NE.ENUM_CENTERED_4TH )  C to be able to exclude this scheme to avoid excessive storage and
565       &   CALL GAD_ADVECTION(bi,bj,tempAdvScheme,GAD_TEMPERATURE,theta,  C recomputation. It *is* differentiable, if you need it.
566       U                      gT,  C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to
567    C disable this section of code.
568            IF (tempMultiDimAdvec) THEN
569    #ifndef DISABLE_DEBUGMODE
570              IF ( debugLevel .GE. debLevB )
571         &     CALL DEBUG_CALL('GAD_ADVECTION',myThid)
572    #endif
573              CALL GAD_ADVECTION(bi,bj,tempAdvScheme,GAD_TEMPERATURE,
574         U                      theta,gT,
575       I                      myTime,myIter,myThid)       I                      myTime,myIter,myThid)
576           IF (saltStepping .AND.          ENDIF
577       &       saltAdvScheme.NE.ENUM_CENTERED_2ND .AND.          IF (saltMultiDimAdvec) THEN
578       &       saltAdvScheme.NE.ENUM_UPWIND_3RD .AND.  #ifndef DISABLE_DEBUGMODE
579       &       saltAdvScheme.NE.ENUM_CENTERED_4TH )            IF ( debugLevel .GE. debLevB )
580       &   CALL GAD_ADVECTION(bi,bj,saltAdvScheme,GAD_SALINITY,salt,       &     CALL DEBUG_CALL('GAD_ADVECTION',myThid)
581       U                      gS,  #endif
582              CALL GAD_ADVECTION(bi,bj,saltAdvScheme,GAD_SALINITY,
583         U                      salt,gS,
584       I                      myTime,myIter,myThid)       I                      myTime,myIter,myThid)
585          ENDIF          ENDIF
586    C Since passive tracers are configurable separately from T,S we
587    C call the multi-dimensional method for PTRACERS regardless
588    C of whether multiDimAdvection is set or not.
589    #ifdef ALLOW_PTRACERS
590            IF ( usePTRACERS ) THEN
591    #ifndef DISABLE_DEBUGMODE
592              IF ( debugLevel .GE. debLevB )
593         &     CALL DEBUG_CALL('PTRACERS_ADVECTION',myThid)
594    #endif
595             CALL PTRACERS_ADVECTION( bi,bj,myIter,myTime,myThid )
596            ENDIF
597    #endif /* ALLOW_PTRACERS */
598    #endif /* DISABLE_MULTIDIM_ADVECTION */
599    
600    #ifndef DISABLE_DEBUGMODE
601           IF ( debugLevel .GE. debLevB )
602         &    CALL DEBUG_MSG('ENTERING DOWNWARD K LOOP',myThid)
603    #endif
604    
605  C--     Start of thermodynamics loop  C--     Start of thermodynamics loop
606          DO k=Nr,1,-1          DO k=Nr,1,-1
# Line 472  C--     Start of thermodynamics loop Line 608  C--     Start of thermodynamics loop
608  C? Patrick Is this formula correct?  C? Patrick Is this formula correct?
609  cph Yes, but I rewrote it.  cph Yes, but I rewrote it.
610  cph Also, the KappaR? need the index and subscript k!  cph Also, the KappaR? need the index and subscript k!
611           kkey = (ikey-1)*Nr + k           kkey = (itdkey-1)*Nr + k
612  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
613    
614  C--       km1    Points to level above k (=k-1)  C--       km1    Points to level above k (=k-1)
# Line 494  C--       Get temporary terms used by te Line 630  C--       Get temporary terms used by te
630       O         xA,yA,uTrans,vTrans,rTrans,maskUp,       O         xA,yA,uTrans,vTrans,rTrans,maskUp,
631       I         myThid)       I         myThid)
632    
633    #ifdef ALLOW_GMREDI
634    
635    C--   Residual transp = Bolus transp + Eulerian transp
636              IF (useGMRedi) THEN
637                CALL GMREDI_CALC_UVFLOW(
638         &            uTrans, vTrans, bi, bj, k, myThid)
639                IF (K.GE.2) CALL GMREDI_CALC_WFLOW(
640         &                    rTrans, bi, bj, k, myThid)
641              ENDIF
642    
643    #ifdef ALLOW_AUTODIFF_TAMC
644    #ifdef GM_BOLUS_ADVEC
645    CADJ STORE uTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
646    CADJ STORE vTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
647    CADJ STORE rTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
648    #endif
649    #endif /* ALLOW_AUTODIFF_TAMC */
650    
651    #endif /* ALLOW_GMREDI */
652    
653  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
654  CADJ STORE KappaRT(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE KappaRT(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte
655  CADJ STORE KappaRS(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE KappaRS(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte
# Line 525  C        and step forward storing result Line 681  C        and step forward storing result
681             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
682       I         bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,       I         bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,
683       I         theta, gT,       I         theta, gT,
      U         gTnm1,  
684       I         myIter, myThid)       I         myIter, myThid)
685           ENDIF           ENDIF
686    
687    #ifdef ALLOW_THERM_SEAICE
688             IF (useThermSeaIce .AND. k.EQ.1) THEN
689               CALL ICE_FREEZE( bi,bj, iMin,iMax,jMin,jMax, myThid )
690             ENDIF
691    #endif
692    
693           IF ( saltStepping ) THEN           IF ( saltStepping ) THEN
694             CALL CALC_GS(             CALL CALC_GS(
695       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
# Line 538  C        and step forward storing result Line 700  C        and step forward storing result
700             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
701       I         bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,       I         bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,
702       I         salt, gS,       I         salt, gS,
      U         gSnm1,  
703       I         myIter, myThid)       I         myIter, myThid)
704           ENDIF           ENDIF
705  #ifdef ALLOW_PASSIVE_TRACER  #ifdef ALLOW_PASSIVE_TRACER
# Line 552  C        and step forward storing result Line 713  C        and step forward storing result
713             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
714       I         bi,bj,iMin,iMax,jMin,jMax,k,tracerAdvScheme,       I         bi,bj,iMin,iMax,jMin,jMax,k,tracerAdvScheme,
715       I         Tr1, gTr1,       I         Tr1, gTr1,
      U         gTr1NM1,  
716       I         myIter,myThid)       I         myIter,myThid)
717           ENDIF           ENDIF
718  #endif  #endif
719    #ifdef ALLOW_PTRACERS
720             IF ( usePTRACERS ) THEN
721               CALL PTRACERS_INTEGRATE(
722         I         bi,bj,k,
723         I         xA,yA,uTrans,vTrans,rTrans,maskUp,
724         X         KappaRS,
725         I         myIter,myTime,myThid)
726             ENDIF
727    #endif /* ALLOW_PTRACERS */
728    
729  #ifdef   ALLOW_OBCS  #ifdef   ALLOW_OBCS
730  C--      Apply open boundary conditions  C--      Apply open boundary conditions
# Line 565  C--      Apply open boundary conditions Line 734  C--      Apply open boundary conditions
734  #endif   /* ALLOW_OBCS */  #endif   /* ALLOW_OBCS */
735    
736  C--      Freeze water  C--      Freeze water
737           IF (allowFreezing) THEN           IF ( allowFreezing .AND. .NOT. useSEAICE
738         &       .AND. .NOT.(useThermSeaIce.AND.k.EQ.1) ) THEN
739  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
740  CADJ STORE gTNm1(:,:,k,bi,bj) = comlev1_bibj_k  CADJ STORE gT(:,:,k,bi,bj) = comlev1_bibj_k
741  CADJ &   , key = kkey, byte = isbyte  CADJ &   , key = kkey, byte = isbyte
742  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
743              CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )              CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )
744           END IF           ENDIF
745    
746  C--     end of thermodynamic k loop (Nr:1)  C--     end of thermodynamic k loop (Nr:1)
747          ENDDO          ENDDO
748    
749    cswdice -- add ---
750    #ifdef ALLOW_THERM_SEAICE
751    c timeaveraging for ice model values
752               CALL ICE_AVE(bi,bj,iMin,iMax,jMin,jMax,myThid )
753    #endif
754    cswdice --- end add ---
755    
756    
757    
 #ifdef ALLOW_AUTODIFF_TAMC  
 C? Patrick? What about this one?  
 cph Keys iikey and idkey dont seem to be needed  
 cph since storing occurs on different tape for each  
 cph impldiff call anyways.  
 cph Thus, common block comlev1_impl isnt needed either.  
 cph Storing below needed in the case useGMREDI.  
         iikey = (ikey-1)*maximpl  
 #endif /* ALLOW_AUTODIFF_TAMC */  
758    
759  C--     Implicit diffusion  C--     Implicit diffusion
760          IF (implicitDiffusion) THEN          IF (implicitDiffusion) THEN
761    
762           IF (tempStepping) THEN           IF (tempStepping) THEN
763  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
764              idkey = iikey + 1  CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
 CADJ STORE gTNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  
765  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
766              CALL IMPLDIFF(              CALL IMPLDIFF(
767       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
# Line 604  CADJ STORE gTNm1(:,:,:,bi,bj) = comlev1_ Line 772  CADJ STORE gTNm1(:,:,:,bi,bj) = comlev1_
772    
773           IF (saltStepping) THEN           IF (saltStepping) THEN
774  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
775           idkey = iikey + 2  CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
 CADJ STORE gSNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  
776  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
777              CALL IMPLDIFF(              CALL IMPLDIFF(
778       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
# Line 617  CADJ STORE gSNm1(:,:,:,bi,bj) = comlev1_ Line 784  CADJ STORE gSNm1(:,:,:,bi,bj) = comlev1_
784  #ifdef ALLOW_PASSIVE_TRACER  #ifdef ALLOW_PASSIVE_TRACER
785           IF (tr1Stepping) THEN           IF (tr1Stepping) THEN
786  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
787  CADJ STORE gTr1Nm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE gTr1(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
788  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
789            CALL IMPLDIFF(            CALL IMPLDIFF(
790       I      bi, bj, iMin, iMax, jMin, jMax,       I      bi, bj, iMin, iMax, jMin, jMax,
# Line 627  CADJ STORE gTr1Nm1(:,:,:,bi,bj) = comlev Line 794  CADJ STORE gTr1Nm1(:,:,:,bi,bj) = comlev
794           ENDIF           ENDIF
795  #endif  #endif
796    
797    #ifdef ALLOW_PTRACERS
798    C Vertical diffusion (implicit) for passive tracers
799             IF ( usePTRACERS ) THEN
800               CALL PTRACERS_IMPLDIFF( bi,bj,KappaRS,myThid )
801             ENDIF
802    #endif /* ALLOW_PTRACERS */
803    
804  #ifdef   ALLOW_OBCS  #ifdef   ALLOW_OBCS
805  C--      Apply open boundary conditions  C--      Apply open boundary conditions
806           IF (useOBCS) THEN           IF (useOBCS) THEN
# Line 639  C--      Apply open boundary conditions Line 813  C--      Apply open boundary conditions
813  C--     End If implicitDiffusion  C--     End If implicitDiffusion
814          ENDIF          ENDIF
815    
816  Ccs-  #ifdef ALLOW_TIMEAVE
817            IF (taveFreq.GT.0. .AND. ivdc_kappa.NE.0.) THEN
818              CALL TIMEAVE_CUMUL_1T(ConvectCountTave, ConvectCount,
819         I                           Nr, deltaTclock, bi, bj, myThid)
820            ENDIF
821            useVariableK = useKPP .OR. useGMredi .OR. ivdc_kappa.NE.0.
822            IF (taveFreq.GT.0. .AND. useVariableK ) THEN
823             IF (implicitDiffusion) THEN
824              CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, gT, kappaRT,
825         I                        Nr, 3, deltaTclock, bi, bj, myThid)
826             ELSE
827              CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, theta, kappaRT,
828         I                        Nr, 3, deltaTclock, bi, bj, myThid)
829             ENDIF
830            ENDIF
831    #endif /* ALLOW_TIMEAVE */
832    
833    #endif /* SINGLE_LAYER_MODE */
834    
835    C--   end bi,bj loops.
836         ENDDO         ENDDO
837        ENDDO        ENDDO
838    
839  #ifdef ALLOW_AIM  #ifndef DISABLE_DEBUGMODE
840        IF ( useAIM ) THEN        If (debugMode) THEN
841         CALL AIM_AIM2DYN_EXCHANGES( myTime, myThid )         CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (THERMODYNAMICS)',myThid)
842           CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (THERMODYNAMICS)',myThid)
843           CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (THERMODYNAMICS)',myThid)
844           CALL DEBUG_STATS_RL(Nr,theta,'Theta (THERMODYNAMICS)',myThid)
845           CALL DEBUG_STATS_RL(Nr,salt,'Salt (THERMODYNAMICS)',myThid)
846           CALL DEBUG_STATS_RL(Nr,Gt,'Gt (THERMODYNAMICS)',myThid)
847           CALL DEBUG_STATS_RL(Nr,Gs,'Gs (THERMODYNAMICS)',myThid)
848           CALL DEBUG_STATS_RL(Nr,GtNm1,'GtNm1 (THERMODYNAMICS)',myThid)
849           CALL DEBUG_STATS_RL(Nr,GsNm1,'GsNm1 (THERMODYNAMICS)',myThid)
850    #ifdef ALLOW_PTRACERS
851           IF ( usePTRACERS ) THEN
852             CALL PTRACERS_DEBUG(myThid)
853           ENDIF
854    #endif /* ALLOW_PTRACERS */
855        ENDIF        ENDIF
856         _EXCH_XYZ_R8(gT,myThid)  #endif
857         _EXCH_XYZ_R8(gS,myThid)  
858  #else  #ifndef DISABLE_DEBUGMODE
859        IF (staggerTimeStep.AND.useCubedSphereExchange) THEN           IF ( debugLevel .GE. debLevB )
860         _EXCH_XYZ_R8(gT,myThid)       &    CALL DEBUG_LEAVE('FORWARD_STEP',myThid)
861         _EXCH_XYZ_R8(gS,myThid)  #endif
       ENDIF  
 #endif /* ALLOW_AIM */  
862    
863        RETURN        RETURN
864        END        END
Legend:
Removed from v.1.10  
changed lines
  Added in v.1.50
  ViewVC Help
Powered by ViewVC 1.1.22