/[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.33 by jmc, Fri Nov 22 03:01:18 2002 UTC revision 1.50 by jmc, Tue Oct 7 04:31:30 2003 UTC
# Line 9  C $Name$ Line 9  C $Name$
9  # ifdef ALLOW_KPP  # ifdef ALLOW_KPP
10  #  include "KPP_OPTIONS.h"  #  include "KPP_OPTIONS.h"
11  # endif  # endif
12  cswdice --- add ----  #ifdef ALLOW_PTRACERS
13  #ifdef ALLOW_THERM_SEAICE  # include "PTRACERS_OPTIONS.h"
 #include "ICE.h"  
14  #endif  #endif
 cswdice ------  
15  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
16    
17  CBOP  CBOP
# Line 84  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"
# Line 96  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 123  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 147  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 159  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    
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
# Line 173  C--   dummy statement to end declaration Line 177  C--   dummy statement to end declaration
177        itdkey = 1        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  
         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 201  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 221  CHPF$&                  ) Line 207  CHPF$&                  )
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 231  C--     Set up work arrays that need val Line 230  C--     Set up work arrays that need val
230            fVerS  (i,j,2) = 0. _d 0            fVerS  (i,j,2) = 0. _d 0
231            fVerTr1(i,j,1) = 0. _d 0            fVerTr1(i,j,1) = 0. _d 0
232            fVerTr1(i,j,2) = 0. _d 0            fVerTr1(i,j,2) = 0. _d 0
           rhoKM1 (i,j)   = 0. _d 0  
233           ENDDO           ENDDO
234          ENDDO          ENDDO
235    
# 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
            phiHyd(i,j,k) = 0. _d 0  
240             sigmaX(i,j,k) = 0. _d 0             sigmaX(i,j,k) = 0. _d 0
241             sigmaY(i,j,k) = 0. _d 0             sigmaY(i,j,k) = 0. _d 0
242             sigmaR(i,j,k) = 0. _d 0             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):
 cph all the following init. are necessary for TAF  
 cph although some of these are re-initialised later.  
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    # 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  # ifdef ALLOW_GMREDI
261             Kwx(i,j,k,bi,bj)  = 0. _d 0             Kwx(i,j,k,bi,bj)  = 0. _d 0
262             Kwy(i,j,k,bi,bj)  = 0. _d 0             Kwy(i,j,k,bi,bj)  = 0. _d 0
# Line 270  cph although some of these are re-initia Line 273  cph although some of these are re-initia
273             GM_PsiX(i,j,k,bi,bj)  = 0. _d 0             GM_PsiX(i,j,k,bi,bj)  = 0. _d 0
274             GM_PsiY(i,j,k,bi,bj)  = 0. _d 0             GM_PsiY(i,j,k,bi,bj)  = 0. _d 0
275  #  endif  #  endif
276    #  ifdef GM_VISBECK_VARIABLE_K
277               VisbeckK(i,j,bi,bj)   = 0. _d 0
278    #  endif
279  # endif /* ALLOW_GMREDI */  # endif /* ALLOW_GMREDI */
280  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
281            ENDDO            ENDDO
# Line 281  cph although some of these are re-initia Line 287  cph although some of these are re-initia
287          jMin = 1-OLy          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=itdkey, byte=isbyte  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
292  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, 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  #ifdef ALLOW_KPP
296  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
297  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
298  #endif  #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 325  C--       Calculate gradients of potenti Line 337  C--       Calculate gradients of potenti
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
 CADJ STORE pressure(:,:,k,bi,bj) =  
 CADJ &     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,       I        bi, bj, iMin, iMax, jMin, jMax, k, k,
# Line 341  CADJ &     comlev1_bibj_k, key=kkey, byt Line 355  CADJ &     comlev1_bibj_k, key=kkey, byt
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
 CADJ STORE pressure(:,:,k-1,bi,bj) =  
 CADJ &     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,       I        bi, bj, iMin, iMax, jMin, jMax, k-1, k,
# Line 350  CADJ &     comlev1_bibj_k, key=kkey, byt Line 362  CADJ &     comlev1_bibj_k, key=kkey, byt
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 364  CADJ STORE rhokm1 (:,:) = comlev1_bibj_k Line 380  CADJ STORE rhokm1 (:,:) = comlev1_bibj_k
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 379  C--     end of diagnostic k loop (Nr:1) Line 399  C--     end of diagnostic k loop (Nr:1)
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=itdkey, byte=isbyte  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
 CADJ STORE pressure (:,:,:,bi,bj) =  
 CADJ &     comlev1_bibj, key=itdkey, byte=isbyte  
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    #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,            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 )
# Line 393  C--     Calculate future values on open Line 415  C--     Calculate future values on open
415  #endif  /* ALLOW_OBCS */  #endif  /* ALLOW_OBCS */
416    
417    
 c********************************************  
 cswdice --- add ---  
418  #ifdef ALLOW_THERM_SEAICE  #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--     including effects from ice  C       including effects from ice
426          CALL ICE_FORCING(          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         ELSE
430    #endif /* ALLOW_THERM_SEAICE */
 cswdice --- end add ---  
431    
432  C--     Determines forcing terms based on external fields  C--     Determines forcing terms based on external fields
433  C       relaxation terms, etc.  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(          CALL EXTERNAL_FORCING_SURF(
439       I             bi, bj, iMin, iMax, jMin, jMax,       I             bi, bj, iMin, iMax, jMin, jMax,
440       I             myThid )       I             myTime, myIter, myThid )
 cswdice --- add ----  
 #endif  
 cswdice --- end add ---  
 c******************************************  
   
   
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
# Line 437  C--     MOST of THERMODYNAMICS will be d Line 463  C--     MOST of THERMODYNAMICS will be d
463  #ifdef  ALLOW_GMREDI  #ifdef  ALLOW_GMREDI
464    
465  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
466  CADJ STORE sigmaX(:,:,k)        = comlev1_bibj_k, key=kkey, byte=isbyte  cph storing here is needed only for one GMREDI_OPTIONS:
467  CADJ STORE sigmaY(:,:,k)        = comlev1_bibj_k, key=kkey, byte=isbyte  cph define GM_BOLUS_ADVEC
468  CADJ STORE sigmaR(:,:,k)        = comlev1_bibj_k, key=kkey, 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    #ifndef DISABLE_DEBUGMODE
477              IF ( debugLevel .GE. debLevB )
478         &     CALL DEBUG_CALL('GMREDI_CALC_TENSOR',myThid)
479    #endif
480            CALL GMREDI_CALC_TENSOR(            CALL GMREDI_CALC_TENSOR(
481       I             bi, bj, iMin, iMax, jMin, jMax,       I             bi, bj, iMin, iMax, jMin, jMax,
482       I             sigmaX, sigmaY, sigmaR,       I             sigmaX, sigmaY, sigmaR,
# Line 467  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 496  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_ Line 534  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_
534  #ifdef ALLOW_PASSIVE_TRACER  #ifdef ALLOW_PASSIVE_TRACER
535  CADJ STORE tr1  (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte  CADJ STORE tr1  (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
536  #endif  #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
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.
546          IF ( useAIM ) THEN          IF ( useAIM ) THEN
547    #ifndef DISABLE_DEBUGMODE
548              IF ( debugLevel .GE. debLevB )
549         &     CALL DEBUG_CALL('AIM_DO_PHYSICS',myThid)
550    #endif
551           CALL TIMER_START('AIM_DO_PHYSICS   [THERMODYNAMICS]', myThid)           CALL TIMER_START('AIM_DO_PHYSICS   [THERMODYNAMICS]', myThid)
552           CALL AIM_DO_PHYSICS( bi, bj, myTime, myIter, myThid )           CALL AIM_DO_PHYSICS( bi, bj, myTime, myIter, myThid )
553           CALL TIMER_STOP( 'AIM_DO_PHYSICS   [THERMODYNAMICS]', myThid)           CALL TIMER_STOP( 'AIM_DO_PHYSICS   [THERMODYNAMICS]', myThid)
554          ENDIF          ENDIF
555  #endif /* ALLOW_AIM */  #endif /* ALLOW_AIM */
556    
 #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 */  
   
557  #ifndef DISABLE_MULTIDIM_ADVECTION  #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.
# Line 526  C recomputation. It *is* differentiable, Line 566  C recomputation. It *is* differentiable,
566  C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to  C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to
567  C disable this section of code.  C disable this section of code.
568          IF (tempMultiDimAdvec) THEN          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,            CALL GAD_ADVECTION(bi,bj,tempAdvScheme,GAD_TEMPERATURE,
574       U                      theta,gT,       U                      theta,gT,
575       I                      myTime,myIter,myThid)       I                      myTime,myIter,myThid)
576          ENDIF          ENDIF
577          IF (saltMultiDimAdvec) THEN          IF (saltMultiDimAdvec) THEN
578    #ifndef DISABLE_DEBUGMODE
579              IF ( debugLevel .GE. debLevB )
580         &     CALL DEBUG_CALL('GAD_ADVECTION',myThid)
581    #endif
582            CALL GAD_ADVECTION(bi,bj,saltAdvScheme,GAD_SALINITY,            CALL GAD_ADVECTION(bi,bj,saltAdvScheme,GAD_SALINITY,
583       U                      salt,gS,       U                      salt,gS,
584       I                      myTime,myIter,myThid)       I                      myTime,myIter,myThid)
# Line 540  C call the multi-dimensional method for Line 588  C call the multi-dimensional method for
588  C of whether multiDimAdvection is set or not.  C of whether multiDimAdvection is set or not.
589  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
590          IF ( usePTRACERS ) THEN          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 )           CALL PTRACERS_ADVECTION( bi,bj,myIter,myTime,myThid )
596          ENDIF          ENDIF
597  #endif /* ALLOW_PTRACERS */  #endif /* ALLOW_PTRACERS */
598  #endif /* DISABLE_MULTIDIM_ADVECTION */  #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
607  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 574  C--       Get temporary terms used by te Line 631  C--       Get temporary terms used by te
631       I         myThid)       I         myThid)
632    
633  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
634    
635  C--   Residual transp = Bolus transp + Eulerian transp  C--   Residual transp = Bolus transp + Eulerian transp
636            IF (useGMRedi) THEN            IF (useGMRedi) THEN
637              CALL GMREDI_CALC_UVFLOW(              CALL GMREDI_CALC_UVFLOW(
# Line 581  C--   Residual transp = Bolus transp + E Line 639  C--   Residual transp = Bolus transp + E
639              IF (K.GE.2) CALL GMREDI_CALC_WFLOW(              IF (K.GE.2) CALL GMREDI_CALC_WFLOW(
640       &                    rTrans, bi, bj, k, myThid)       &                    rTrans, bi, bj, k, myThid)
641            ENDIF            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 */  #endif /* ALLOW_GMREDI */
652    
653  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 616  C        and step forward storing result Line 683  C        and step forward storing result
683       I         theta, gT,       I         theta, gT,
684       I         myIter, myThid)       I         myIter, myThid)
685           ENDIF           ENDIF
686  cswdice ---- add ---  
687  #ifdef ALLOW_THERM_SEAICE  #ifdef ALLOW_THERM_SEAICE
688         if (k.eq.1) then           IF (useThermSeaIce .AND. k.EQ.1) THEN
689          call ICE_FREEZE(bi, bj, iMin, iMax, jMin, jMax, myThid )             CALL ICE_FREEZE( bi,bj, iMin,iMax,jMin,jMax, myThid )
690         endif           ENDIF
691  #endif  #endif
692  cswdice -- end add ---  
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 651  cswdice -- end add --- Line 718  cswdice -- end add ---
718  #endif  #endif
719  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
720           IF ( usePTRACERS ) THEN           IF ( usePTRACERS ) THEN
721             CALL PTRACERS_INTEGERATE(             CALL PTRACERS_INTEGRATE(
722       I         bi,bj,k,       I         bi,bj,k,
723       I         xA,yA,uTrans,vTrans,rTrans,maskUp,       I         xA,yA,uTrans,vTrans,rTrans,maskUp,
724       X         KappaRS,       X         KappaRS,
# Line 667  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 .AND. .NOT. useSEAICE ) 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 gT(:,:,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
# Line 745  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    #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 */  #endif /* SINGLE_LAYER_MODE */
834    
835  Ccs-  C--   end bi,bj loops.
836         ENDDO         ENDDO
837        ENDDO        ENDDO
838    
 #ifdef ALLOW_AIM  
 c     IF ( useAIM ) THEN  
 c      CALL AIM_AIM2DYN_EXCHANGES( myTime, myThid )  
 c     ENDIF  
 #endif /* ALLOW_AIM */  
 c     IF ( staggerTimeStep ) THEN  
 c      IF ( useAIM .OR. useCubedSphereExchange ) THEN  
 c        IF (tempStepping) _EXCH_XYZ_R8(gT,myThid)  
 c        IF (saltStepping) _EXCH_XYZ_R8(gS,myThid)  
 c      ELSEIF ( useGMRedi .AND. Oly.LT.4 ) THEN  
 cc       .AND. GM_AdvForm .AND. .NOT.GM_AdvSeparate ) THEN  
 c        IF (tempMultiDimAdvec) _EXCH_XYZ_R8(gT,myThid)  
 c        IF (saltMultiDimAdvec) _EXCH_XYZ_R8(gS,myThid)  
 c      ENDIF    
 c     ENDIF    
   
839  #ifndef DISABLE_DEBUGMODE  #ifndef DISABLE_DEBUGMODE
840        If (debugMode) THEN        If (debugMode) THEN
841         CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (THERMODYNAMICS)',myThid)         CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (THERMODYNAMICS)',myThid)
# Line 786  c     ENDIF Line 855  c     ENDIF
855        ENDIF        ENDIF
856  #endif  #endif
857    
858    #ifndef DISABLE_DEBUGMODE
859             IF ( debugLevel .GE. debLevB )
860         &    CALL DEBUG_LEAVE('FORWARD_STEP',myThid)
861    #endif
862    
863        RETURN        RETURN
864        END        END
Legend:
Removed from v.1.33  
changed lines
  Added in v.1.50
  ViewVC Help
Powered by ViewVC 1.1.22