/[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.5 by heimbach, Mon Sep 10 16:35:27 2001 UTC revision 1.25 by heimbach, Sat Jul 13 04:59:42 2002 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    #endif /* ALLOW_AUTODIFF_TAMC */
13    
14    CBOP
15    C     !ROUTINE: THERMODYNAMICS
16    C     !INTERFACE:
17        SUBROUTINE THERMODYNAMICS(myTime, myIter, myThid)        SUBROUTINE THERMODYNAMICS(myTime, myIter, myThid)
18  C     /==========================================================\  C     !DESCRIPTION: \bv
19  C     | SUBROUTINE THERMODYNAMICS                                |  C     *==========================================================*
20  C     | o Controlling routine for the prognostic part of the     |  C     | SUBROUTINE THERMODYNAMICS                                
21  C     |   thermo-dynamics.                                       |  C     | o Controlling routine for the prognostic part of the      
22  C     |==========================================================|  C     |   thermo-dynamics.                                        
23  C     \==========================================================/  C     *===========================================================
24        IMPLICIT NONE  C     | The algorithm...
25    C     |
26    C     | "Correction Step"
27    C     | =================
28    C     | Here we update the horizontal velocities with the surface
29    C     | pressure such that the resulting flow is either consistent
30    C     | with the free-surface evolution or the rigid-lid:
31    C     |   U[n] = U* + dt x d/dx P
32    C     |   V[n] = V* + dt x d/dy P
33    C     |
34    C     | "Calculation of Gs"
35    C     | ===================
36    C     | This is where all the accelerations and tendencies (ie.
37    C     | physics, parameterizations etc...) are calculated
38    C     |   rho = rho ( theta[n], salt[n] )
39    C     |   b   = b(rho, theta)
40    C     |   K31 = K31 ( rho )
41    C     |   Gu[n] = Gu( u[n], v[n], wVel, b, ... )
42    C     |   Gv[n] = Gv( u[n], v[n], wVel, b, ... )
43    C     |   Gt[n] = Gt( theta[n], u[n], v[n], wVel, K31, ... )
44    C     |   Gs[n] = Gs( salt[n], u[n], v[n], wVel, K31, ... )
45    C     |
46    C     | "Time-stepping" or "Prediction"
47    C     | ================================
48    C     | The models variables are stepped forward with the appropriate
49    C     | time-stepping scheme (currently we use Adams-Bashforth II)
50    C     | - For momentum, the result is always *only* a "prediction"
51    C     | in that the flow may be divergent and will be "corrected"
52    C     | later with a surface pressure gradient.
53    C     | - Normally for tracers the result is the new field at time
54    C     | level [n+1} *BUT* in the case of implicit diffusion the result
55    C     | is also *only* a prediction.
56    C     | - We denote "predictors" with an asterisk (*).
57    C     |   U* = U[n] + dt x ( 3/2 Gu[n] - 1/2 Gu[n-1] )
58    C     |   V* = V[n] + dt x ( 3/2 Gv[n] - 1/2 Gv[n-1] )
59    C     |   theta[n+1] = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
60    C     |   salt[n+1] = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
61    C     | With implicit diffusion:
62    C     |   theta* = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
63    C     |   salt* = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
64    C     |   (1 + dt * K * d_zz) theta[n] = theta*
65    C     |   (1 + dt * K * d_zz) salt[n] = salt*
66    C     |
67    C     *==========================================================*
68    C     \ev
69    
70    C     !USES:
71          IMPLICIT NONE
72  C     == Global variables ===  C     == Global variables ===
73  #include "SIZE.h"  #include "SIZE.h"
74  #include "EEPARAMS.h"  #include "EEPARAMS.h"
# Line 22  C     == Global variables === Line 79  C     == Global variables ===
79  #ifdef ALLOW_PASSIVE_TRACER  #ifdef ALLOW_PASSIVE_TRACER
80  #include "TR1.h"  #include "TR1.h"
81  #endif  #endif
   
82  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
83  # include "tamc.h"  # include "tamc.h"
84  # include "tamc_keys.h"  # include "tamc_keys.h"
# Line 34  C     == Global variables === Line 90  C     == Global variables ===
90  #  include "GMREDI.h"  #  include "GMREDI.h"
91  # endif  # endif
92  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
   
93  #ifdef ALLOW_TIMEAVE  #ifdef ALLOW_TIMEAVE
94  #include "TIMEAVE_STATV.h"  #include "TIMEAVE_STATV.h"
95  #endif  #endif
96    
97    C     !INPUT/OUTPUT PARAMETERS:
98  C     == Routine arguments ==  C     == Routine arguments ==
99  C     myTime - Current time in simulation  C     myTime - Current time in simulation
100  C     myIter - Current iteration number in simulation  C     myIter - Current iteration number in simulation
# Line 47  C     myThid - Thread number for this in Line 103  C     myThid - Thread number for this in
103        INTEGER myIter        INTEGER myIter
104        INTEGER myThid        INTEGER myThid
105    
106    C     !LOCAL VARIABLES:
107  C     == Local variables  C     == Local variables
108  C     xA, yA                 - Per block temporaries holding face areas  C     xA, yA                 - Per block temporaries holding face areas
109  C     uTrans, vTrans, rTrans - Per block temporaries holding flow  C     uTrans, vTrans, rTrans - Per block temporaries holding flow
# Line 94  C                      index into fVerTe Line 151  C                      index into fVerTe
151        _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
152        _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
153        _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
154    C     This is currently used by IVDC and Diagnostics
 C This is currently used by IVDC and Diagnostics  
155        _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
   
156        INTEGER iMin, iMax        INTEGER iMin, iMax
157        INTEGER jMin, jMax        INTEGER jMin, jMax
158        INTEGER bi, bj        INTEGER bi, bj
159        INTEGER i, j        INTEGER i, j
160        INTEGER k, km1, kup, kDown        INTEGER k, km1, kup, kDown
161    
162  Cjmc : add for phiHyd output <- but not working if multi tile per CPU  CEOP
 c     CHARACTER*(MAX_LEN_MBUF) suff  
 c     LOGICAL  DIFFERENT_MULTIPLE  
 c     EXTERNAL DIFFERENT_MULTIPLE  
 Cjmc(end)  
163    
 C---    The algorithm...  
 C  
 C       "Correction Step"  
 C       =================  
 C       Here we update the horizontal velocities with the surface  
 C       pressure such that the resulting flow is either consistent  
 C       with the free-surface evolution or the rigid-lid:  
 C         U[n] = U* + dt x d/dx P  
 C         V[n] = V* + dt x d/dy P  
 C  
 C       "Calculation of Gs"  
 C       ===================  
 C       This is where all the accelerations and tendencies (ie.  
 C       physics, parameterizations etc...) are calculated  
 C         rho = rho ( theta[n], salt[n] )  
 C         b   = b(rho, theta)  
 C         K31 = K31 ( rho )  
 C         Gu[n] = Gu( u[n], v[n], wVel, b, ... )  
 C         Gv[n] = Gv( u[n], v[n], wVel, b, ... )  
 C         Gt[n] = Gt( theta[n], u[n], v[n], wVel, K31, ... )  
 C         Gs[n] = Gs( salt[n], u[n], v[n], wVel, K31, ... )  
 C  
 C       "Time-stepping" or "Prediction"  
 C       ================================  
 C       The models variables are stepped forward with the appropriate  
 C       time-stepping scheme (currently we use Adams-Bashforth II)  
 C       - For momentum, the result is always *only* a "prediction"  
 C       in that the flow may be divergent and will be "corrected"  
 C       later with a surface pressure gradient.  
 C       - Normally for tracers the result is the new field at time  
 C       level [n+1} *BUT* in the case of implicit diffusion the result  
 C       is also *only* a prediction.  
 C       - We denote "predictors" with an asterisk (*).  
 C         U* = U[n] + dt x ( 3/2 Gu[n] - 1/2 Gu[n-1] )  
 C         V* = V[n] + dt x ( 3/2 Gv[n] - 1/2 Gv[n-1] )  
 C         theta[n+1] = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )  
 C         salt[n+1] = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )  
 C       With implicit diffusion:  
 C         theta* = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )  
 C         salt* = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )  
 C         (1 + dt * K * d_zz) theta[n] = theta*  
 C         (1 + dt * K * d_zz) salt[n] = salt*  
 C---  
   
164  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
165  C--   dummy statement to end declaration part  C--   dummy statement to end declaration part
166        ikey = 1        ikey = 1
# Line 170  C     uninitialised but inert locations. Line 177  C     uninitialised but inert locations.
177          yA(i,j)      = 0. _d 0          yA(i,j)      = 0. _d 0
178          uTrans(i,j)  = 0. _d 0          uTrans(i,j)  = 0. _d 0
179          vTrans(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  
180          rhok   (i,j) = 0. _d 0          rhok   (i,j) = 0. _d 0
181          phiSurfX(i,j) = 0. _d 0          phiSurfX(i,j) = 0. _d 0
182          phiSurfY(i,j) = 0. _d 0          phiSurfY(i,j) = 0. _d 0
# Line 204  CHPF$&                  ) Line 204  CHPF$&                  )
204  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
205            act1 = bi - myBxLo(myThid)            act1 = bi - myBxLo(myThid)
206            max1 = myBxHi(myThid) - myBxLo(myThid) + 1            max1 = myBxHi(myThid) - myBxLo(myThid) + 1
   
207            act2 = bj - myByLo(myThid)            act2 = bj - myByLo(myThid)
208            max2 = myByHi(myThid) - myByLo(myThid) + 1            max2 = myByHi(myThid) - myByLo(myThid) + 1
   
209            act3 = myThid - 1            act3 = myThid - 1
210            max3 = nTx*nTy            max3 = nTx*nTy
   
211            act4 = ikey_dynamics - 1            act4 = ikey_dynamics - 1
   
212            ikey = (act1 + 1) + act2*max1            ikey = (act1 + 1) + act2*max1
213       &                      + act3*max1*max2       &                      + act3*max1*max2
214       &                      + act4*max1*max2*max3       &                      + act4*max1*max2*max3
# Line 228  C--     Set up work arrays that need val Line 224  C--     Set up work arrays that need val
224            fVerS  (i,j,2) = 0. _d 0            fVerS  (i,j,2) = 0. _d 0
225            fVerTr1(i,j,1) = 0. _d 0            fVerTr1(i,j,1) = 0. _d 0
226            fVerTr1(i,j,2) = 0. _d 0            fVerTr1(i,j,2) = 0. _d 0
227              rhoKM1 (i,j)   = 0. _d 0
228           ENDDO           ENDDO
229          ENDDO          ENDDO
230    
# Line 235  C--     Set up work arrays that need val Line 232  C--     Set up work arrays that need val
232           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
233            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
234  C This is currently also used by IVDC and Diagnostics  C This is currently also used by IVDC and Diagnostics
235               phiHyd(i,j,k)  = 0. _d 0
236               sigmaX(i,j,k) = 0. _d 0
237               sigmaY(i,j,k) = 0. _d 0
238               sigmaR(i,j,k) = 0. _d 0
239             ConvectCount(i,j,k) = 0.             ConvectCount(i,j,k) = 0.
240             KappaRT(i,j,k) = 0. _d 0             KappaRT(i,j,k) = 0. _d 0
241             KappaRS(i,j,k) = 0. _d 0             KappaRS(i,j,k) = 0. _d 0
# Line 244  C This is currently also used by IVDC an Line 245  C This is currently also used by IVDC an
245  #ifdef ALLOW_PASSIVE_TRACER  #ifdef ALLOW_PASSIVE_TRACER
246             gTr1(i,j,k,bi,bj) = 0. _d 0             gTr1(i,j,k,bi,bj) = 0. _d 0
247  #endif  #endif
248    #ifdef ALLOW_GMREDI
249               Kwx(i,j,k,bi,bj)    = 0. _d 0
250               Kwy(i,j,k,bi,bj)    = 0. _d 0
251               Kwz(i,j,k,bi,bj)    = 0. _d 0
252    #ifdef GM_NON_UNITY_DIAGONAL
253               Kux(i,j,k,bi,bj)    = 0. _d 0
254               Kvy(i,j,k,bi,bj)    = 0. _d 0
255    #endif
256    #endif /* ALLOW_GMREDI */
257  #endif  #endif
258            ENDDO            ENDDO
259           ENDDO           ENDDO
260          ENDDO          ENDDO
261    
262          iMin = 1-OLx+1          iMin = 1-OLx
263          iMax = sNx+OLx          iMax = sNx+OLx
264          jMin = 1-OLy+1          jMin = 1-OLy
265          jMax = sNy+OLy          jMax = sNy+OLy
266    
267    
268  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
269  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
270  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
271    #ifdef ALLOW_KPP
272    CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
273    CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
274    #endif
275  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
276    
277  C--     Start of diagnostic loop  C--     Start of diagnostic loop
# Line 288  C--       Apply OBC to W if in N-H mode Line 302  C--       Apply OBC to W if in N-H mode
302  #endif    /* ALLOW_NONHYDROSTATIC */  #endif    /* ALLOW_NONHYDROSTATIC */
303  #endif    /* ALLOW_OBCS */  #endif    /* ALLOW_OBCS */
304    
305    C--       Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h
306    C--       MOST of THERMODYNAMICS will be disabled
307    #ifndef SINGLE_LAYER_MODE
308    
309  C--       Calculate gradients of potential density for isoneutral  C--       Calculate gradients of potential density for isoneutral
310  C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)  C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)
311  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
# Line 329  c ==> should use sigmaR !!! Line 347  c ==> should use sigmaR !!!
347       I        myTime, myIter, myThid)       I        myTime, myIter, myThid)
348            ENDIF            ENDIF
349    
350    #endif /* SINGLE_LAYER_MODE */
351    
352  C--     end of diagnostic k loop (Nr:1)  C--     end of diagnostic k loop (Nr:1)
353          ENDDO          ENDDO
354    
355  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
356  cph avoids recomputation of integrate_for_w  cph avoids recomputation of integrate_for_w
357  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
358  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
359    
360  #ifdef  ALLOW_OBCS  #ifdef  ALLOW_OBCS
361  C--     Calculate future values on open boundaries  C--     Calculate future values on open boundaries
362          IF (useOBCS) THEN          IF (useOBCS) THEN
363            CALL OBCS_CALC( bi, bj, myTime+deltaT,            CALL OBCS_CALC( bi, bj, myTime+deltaT, myIter+1,
364       I            uVel, vVel, wVel, theta, salt,       I            uVel, vVel, wVel, theta, salt,
365       I            myThid )       I            myThid )
366          ENDIF          ENDIF
# Line 363  CADJ STORE surfacetendencyT(:,:,bi,bj) Line 383  CADJ STORE surfacetendencyT(:,:,bi,bj)
383  CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  CADJ &     = comlev1_bibj, key=ikey, byte=isbyte
384  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
385    
386    C--     Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h
387    C--     MOST of THERMODYNAMICS will be disabled
388    #ifndef SINGLE_LAYER_MODE
389    
390  #ifdef  ALLOW_GMREDI  #ifdef  ALLOW_GMREDI
391    
392  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
393  CADJ STORE sigmaX(:,:,:) = comlev1, key=ikey, byte=isbyte  CADJ STORE sigmaX(:,:,k)        = comlev1_bibj_k, key=kkey, byte=isbyte
394  CADJ STORE sigmaY(:,:,:) = comlev1, key=ikey, byte=isbyte  CADJ STORE sigmaY(:,:,k)        = comlev1_bibj_k, key=kkey, byte=isbyte
395  CADJ STORE sigmaR(:,:,:) = comlev1, key=ikey, byte=isbyte  CADJ STORE sigmaR(:,:,k)        = comlev1_bibj_k, key=kkey, byte=isbyte
396  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
397  C--     Calculate iso-neutral slopes for the GM/Redi parameterisation  C--     Calculate iso-neutral slopes for the GM/Redi parameterisation
398          IF (useGMRedi) THEN          IF (useGMRedi) THEN
399            DO k=1,Nr            CALL GMREDI_CALC_TENSOR(
400              CALL GMREDI_CALC_TENSOR(       I             bi, bj, iMin, iMax, jMin, jMax,
      I             bi, bj, iMin, iMax, jMin, jMax, k,  
401       I             sigmaX, sigmaY, sigmaR,       I             sigmaX, sigmaY, sigmaR,
402       I             myThid )       I             myThid )
           ENDDO  
403  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
404          ELSE          ELSE
405            DO k=1, Nr            CALL GMREDI_CALC_TENSOR_DUMMY(
406              CALL GMREDI_CALC_TENSOR_DUMMY(       I             bi, bj, iMin, iMax, jMin, jMax,
      I             bi, bj, iMin, iMax, jMin, jMax, k,  
407       I             sigmaX, sigmaY, sigmaR,       I             sigmaX, sigmaY, sigmaR,
408       I             myThid )       I             myThid )
           ENDDO  
409  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
410          ENDIF          ENDIF
411    
# Line 411  C--     Compute KPP mixing coefficients Line 431  C--     Compute KPP mixing coefficients
431    
432  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
433  CADJ STORE KPPghat   (:,:,:,bi,bj)  CADJ STORE KPPghat   (:,:,:,bi,bj)
 CADJ &   , KPPviscAz (:,:,:,bi,bj)  
434  CADJ &   , KPPdiffKzT(:,:,:,bi,bj)  CADJ &   , KPPdiffKzT(:,:,:,bi,bj)
435  CADJ &   , KPPdiffKzS(:,:,:,bi,bj)  CADJ &   , KPPdiffKzS(:,:,:,bi,bj)
436  CADJ &   , KPPfrac   (:,:  ,bi,bj)  CADJ &   , KPPfrac   (:,:  ,bi,bj)
# Line 421  CADJ &                 = comlev1_bibj, k Line 440  CADJ &                 = comlev1_bibj, k
440  #endif  /* ALLOW_KPP */  #endif  /* ALLOW_KPP */
441    
442  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
443  CADJ STORE KappaRT(:,:,:)     = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE KappaRT(:,:,:)     = comlev1_bibj, key=ikey, byte=isbyte
444  CADJ STORE KappaRS(:,:,:)     = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE KappaRS(:,:,:)     = comlev1_bibj, key=ikey, byte=isbyte
445  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
446  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
447  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
448  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
449  #ifdef ALLOW_PASSIVE_TRACER  #ifdef ALLOW_PASSIVE_TRACER
450  CADJ STORE tr1  (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE tr1  (:,:,:,bi,bj) = comlev1_bibj, key=ikey, byte=isbyte
451  #endif  #endif
452  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
453    
# Line 442  C note(jmc) : phiHyd=0 at this point but Line 461  C note(jmc) : phiHyd=0 at this point but
461          ENDIF          ENDIF
462  #endif /* ALLOW_AIM */  #endif /* ALLOW_AIM */
463    
464    #ifdef ALLOW_TIMEAVE
465            IF (taveFreq.GT.0. .AND. ivdc_kappa.NE.0.) THEN
466              CALL TIMEAVE_CUMULATE(ConvectCountTave, ConvectCount, Nr,
467         I                               deltaTclock, bi, bj, myThid)
468            ENDIF
469    #endif /* ALLOW_TIMEAVE */
470    
471    #ifndef DISABLE_MULTIDIM_ADVECTION
472  C--     Some advection schemes are better calculated using a multi-dimensional  C--     Some advection schemes are better calculated using a multi-dimensional
473  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.
474  cph        IF (multiDimAdvection) THEN  C
475  cph         IF (tempStepping .AND.  C The CPP flag DISABLE_MULTIDIM_ADVECTION is currently unset in GAD_OPTIONS.h
476  cph     &       tempAdvScheme.NE.ENUM_CENTERED_2ND .AND.  C The default is to use multi-dimensinal advection for non-linear advection
477  cph     &       tempAdvScheme.NE.ENUM_UPWIND_3RD .AND.  C schemes. However, for the sake of efficiency of the adjoint it is necessary
478  cph     &       tempAdvScheme.NE.ENUM_CENTERED_4TH )  C to be able to exclude this scheme to avoid excessive storage and
479  cph     &   CALL GAD_ADVECTION(bi,bj,tempAdvScheme,GAD_TEMPERATURE,theta,  C recomputation. It *is* differentiable, if you need it.
480  cph     U                      gT,  C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to
481  cph     I                      myTime,myIter,myThid)  C disable this section of code.
482  cph         IF (saltStepping .AND.          IF (tempMultiDimAdvec) THEN
483  cph     &       saltAdvScheme.NE.ENUM_CENTERED_2ND .AND.            CALL GAD_ADVECTION(bi,bj,tempAdvScheme,GAD_TEMPERATURE,
484  cph     &       saltAdvScheme.NE.ENUM_UPWIND_3RD .AND.       U                      theta,gT,
485  cph     &       saltAdvScheme.NE.ENUM_CENTERED_4TH )       I                      myTime,myIter,myThid)
486  cph     &   CALL GAD_ADVECTION(bi,bj,saltAdvScheme,GAD_SALINITY,salt,          ENDIF
487  cph     U                      gS,          IF (saltMultiDimAdvec) THEN
488  cph     I                      myTime,myIter,myThid)            CALL GAD_ADVECTION(bi,bj,saltAdvScheme,GAD_SALINITY,
489  cph        ENDIF       U                      salt,gS,
490         I                      myTime,myIter,myThid)
491            ENDIF
492    C Since passive tracers are configurable separately from T,S we
493    C call the multi-dimensional method for PTRACERS regardless
494    C of whether multiDimAdvection is set or not.
495    #ifdef ALLOW_PTRACERS
496            IF ( usePTRACERS ) THEN
497             CALL PTRACERS_ADVECTION( bi,bj,myIter,myTime,myThid )
498            ENDIF
499    #endif /* ALLOW_PTRACERS */
500    #endif /* DISABLE_MULTIDIM_ADVECTION */
501    
502  C--     Start of thermodynamics loop  C--     Start of thermodynamics loop
503          DO k=Nr,1,-1          DO k=Nr,1,-1
# Line 490  C--       Get temporary terms used by te Line 527  C--       Get temporary terms used by te
527       O         xA,yA,uTrans,vTrans,rTrans,maskUp,       O         xA,yA,uTrans,vTrans,rTrans,maskUp,
528       I         myThid)       I         myThid)
529    
530    #ifdef ALLOW_GMREDI
531    C--   Residual transp = Bolus transp + Eulerian transp
532              IF (useGMRedi) THEN
533                CALL GMREDI_CALC_UVFLOW(
534         &            uTrans, vTrans, bi, bj, k, myThid)
535                IF (K.GE.2) CALL GMREDI_CALC_WFLOW(
536         &                    rTrans, bi, bj, k, myThid)
537              ENDIF
538    #endif /* ALLOW_GMREDI */
539    
540  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
541  CADJ STORE KappaRT(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE KappaRT(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte
542  CADJ STORE KappaRS(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE KappaRS(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte
# Line 517  C        and step forward storing result Line 564  C        and step forward storing result
564       I         xA,yA,uTrans,vTrans,rTrans,maskUp,       I         xA,yA,uTrans,vTrans,rTrans,maskUp,
565       I         KappaRT,       I         KappaRT,
566       U         fVerT,       U         fVerT,
567       I         myTime, myThid)       I         myTime,myIter,myThid)
568             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
569       I         bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,       I         bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,
570       I         theta, gT,       I         theta, gT,
      U         gTnm1,  
571       I         myIter, myThid)       I         myIter, myThid)
572           ENDIF           ENDIF
573           IF ( saltStepping ) THEN           IF ( saltStepping ) THEN
# Line 530  C        and step forward storing result Line 576  C        and step forward storing result
576       I         xA,yA,uTrans,vTrans,rTrans,maskUp,       I         xA,yA,uTrans,vTrans,rTrans,maskUp,
577       I         KappaRS,       I         KappaRS,
578       U         fVerS,       U         fVerS,
579       I         myTime, myThid)       I         myTime,myIter,myThid)
580             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
581       I         bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,       I         bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,
582       I         salt, gS,       I         salt, gS,
      U         gSnm1,  
583       I         myIter, myThid)       I         myIter, myThid)
584           ENDIF           ENDIF
585  #ifdef ALLOW_PASSIVE_TRACER  #ifdef ALLOW_PASSIVE_TRACER
# Line 544  C        and step forward storing result Line 589  C        and step forward storing result
589       I         xA,yA,uTrans,vTrans,rTrans,maskUp,       I         xA,yA,uTrans,vTrans,rTrans,maskUp,
590       I         KappaRT,       I         KappaRT,
591       U         fVerTr1,       U         fVerTr1,
592       I         myTime, myThid)       I         myTime,myIter,myThid)
593             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
594       I         bi,bj,iMin,iMax,jMin,jMax,k,tracerAdvScheme,       I         bi,bj,iMin,iMax,jMin,jMax,k,tracerAdvScheme,
595       I         Tr1, gTr1,       I         Tr1, gTr1,
596       U         gTr1NM1,       I         myIter,myThid)
      I         myIter, myThid)  
597           ENDIF           ENDIF
598  #endif  #endif
599    #ifdef ALLOW_PTRACERS
600             IF ( usePTRACERS ) THEN
601               CALL PTRACERS_INTEGERATE(
602         I         bi,bj,k,
603         I         xA,yA,uTrans,vTrans,rTrans,maskUp,
604         X         KappaRS,
605         I         myIter,myTime,myThid)
606             ENDIF
607    #endif /* ALLOW_PTRACERS */
608    
609  #ifdef   ALLOW_OBCS  #ifdef   ALLOW_OBCS
610  C--      Apply open boundary conditions  C--      Apply open boundary conditions
611           IF (useOBCS) THEN           IF (useOBCS) THEN
612             CALL OBCS_APPLY_TS( bi, bj, k, gTnm1, gSnm1, myThid )             CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )
613           END IF           END IF
614  #endif   /* ALLOW_OBCS */  #endif   /* ALLOW_OBCS */
615    
616  C--      Freeze water  C--      Freeze water
617           IF (allowFreezing) THEN           IF (allowFreezing) THEN
618  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
619  CADJ STORE gTNm1(:,:,k,bi,bj) = comlev1_bibj_k  CADJ STORE gT(:,:,k,bi,bj) = comlev1_bibj_k
620  CADJ &   , key = kkey, byte = isbyte  CADJ &   , key = kkey, byte = isbyte
621  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
622              CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )              CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )
# Line 572  CADJ &   , key = kkey, byte = isbyte Line 625  CADJ &   , key = kkey, byte = isbyte
625  C--     end of thermodynamic k loop (Nr:1)  C--     end of thermodynamic k loop (Nr:1)
626          ENDDO          ENDDO
627    
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 C? Patrick? What about this one?  
 cph Keys iikey and idkey don't seem to be needed  
 cph since storing occurs on different tape for each  
 cph impldiff call anyways.  
 cph Thus, common block comlev1_impl isn't needed either.  
 cph Storing below needed in the case useGMREDI.  
         iikey = (ikey-1)*maximpl  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
628  C--     Implicit diffusion  C--     Implicit diffusion
629          IF (implicitDiffusion) THEN          IF (implicitDiffusion) THEN
630    
631           IF (tempStepping) THEN           IF (tempStepping) THEN
632  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
633              idkey = iikey + 1  CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
 CADJ STORE gTNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  
634  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
635              CALL IMPLDIFF(              CALL IMPLDIFF(
636       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
637       I         deltaTtracer, KappaRT, recip_HFacC,       I         deltaTtracer, KappaRT, recip_HFacC,
638       U         gTNm1,       U         gT,
639       I         myThid )       I         myThid )
640           ENDIF           ENDIF
641    
642           IF (saltStepping) THEN           IF (saltStepping) THEN
643  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
644           idkey = iikey + 2  CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
 CADJ STORE gSNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  
645  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
646              CALL IMPLDIFF(              CALL IMPLDIFF(
647       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
648       I         deltaTtracer, KappaRS, recip_HFacC,       I         deltaTtracer, KappaRS, recip_HFacC,
649       U         gSNm1,       U         gS,
650       I         myThid )       I         myThid )
651           ENDIF           ENDIF
652    
653  #ifdef ALLOW_PASSIVE_TRACER  #ifdef ALLOW_PASSIVE_TRACER
654           IF (tr1Stepping) THEN           IF (tr1Stepping) THEN
655  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
656  CADJ STORE gTr1Nm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE gTr1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte
657  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
658            CALL IMPLDIFF(            CALL IMPLDIFF(
659       I      bi, bj, iMin, iMax, jMin, jMax,       I      bi, bj, iMin, iMax, jMin, jMax,
660       I      deltaTtracer, KappaRT, recip_HFacC,       I      deltaTtracer, KappaRT, recip_HFacC,
661       U      gTr1Nm1,       U      gTr1,
662       I      myThid )       I      myThid )
663           ENDIF           ENDIF
664  #endif  #endif
665    
666    #ifdef ALLOW_PTRACERS
667    C Vertical diffusion (implicit) for passive tracers
668             IF ( usePTRACERS ) THEN
669               CALL PTRACERS_IMPLDIFF( bi,bj,KappaRS,myThid )
670             ENDIF
671    #endif /* ALLOW_PTRACERS */
672    
673  #ifdef   ALLOW_OBCS  #ifdef   ALLOW_OBCS
674  C--      Apply open boundary conditions  C--      Apply open boundary conditions
675           IF (useOBCS) THEN           IF (useOBCS) THEN
676             DO K=1,Nr             DO K=1,Nr
677               CALL OBCS_APPLY_TS( bi, bj, k, gTnm1, gSnm1, myThid )               CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )
678             ENDDO             ENDDO
679           END IF           END IF
680  #endif   /* ALLOW_OBCS */  #endif   /* ALLOW_OBCS */
# Line 635  C--      Apply open boundary conditions Line 682  C--      Apply open boundary conditions
682  C--     End If implicitDiffusion  C--     End If implicitDiffusion
683          ENDIF          ENDIF
684    
685    #endif /* SINGLE_LAYER_MODE */
686    
687  Ccs-  Ccs-
688         ENDDO         ENDDO
689        ENDDO        ENDDO
# Line 643  Ccs- Line 692  Ccs-
692        IF ( useAIM ) THEN        IF ( useAIM ) THEN
693         CALL AIM_AIM2DYN_EXCHANGES( myTime, myThid )         CALL AIM_AIM2DYN_EXCHANGES( myTime, myThid )
694        ENDIF        ENDIF
        _EXCH_XYZ_R8(gTnm1,myThid)  
        _EXCH_XYZ_R8(gSnm1,myThid)  
 #else  
       IF (staggerTimeStep.AND.useCubedSphereExchange) THEN  
        _EXCH_XYZ_R8(gTnm1,myThid)  
        _EXCH_XYZ_R8(gSnm1,myThid)  
       ENDIF  
695  #endif /* ALLOW_AIM */  #endif /* ALLOW_AIM */
696          IF ( staggerTimeStep ) THEN
697           IF ( useAIM .OR. useCubedSphereExchange ) THEN
698             IF (tempStepping) _EXCH_XYZ_R8(gT,myThid)
699             IF (saltStepping) _EXCH_XYZ_R8(gS,myThid)
700           ELSEIF ( useGMRedi .AND. Oly.LT.4 ) THEN
701    c        .AND. GM_AdvForm .AND. .NOT.GM_AdvSeparate ) THEN
702             IF (tempMultiDimAdvec) _EXCH_XYZ_R8(gT,myThid)
703             IF (saltMultiDimAdvec) _EXCH_XYZ_R8(gS,myThid)
704           ENDIF  
705          ENDIF  
706    
707    #ifndef DISABLE_DEBUGMODE
708          If (debugMode) THEN
709           CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (THERMODYNAMICS)',myThid)
710           CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (THERMODYNAMICS)',myThid)
711           CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (THERMODYNAMICS)',myThid)
712           CALL DEBUG_STATS_RL(Nr,theta,'Theta (THERMODYNAMICS)',myThid)
713           CALL DEBUG_STATS_RL(Nr,salt,'Salt (THERMODYNAMICS)',myThid)
714           CALL DEBUG_STATS_RL(Nr,Gt,'Gt (THERMODYNAMICS)',myThid)
715           CALL DEBUG_STATS_RL(Nr,Gs,'Gs (THERMODYNAMICS)',myThid)
716           CALL DEBUG_STATS_RL(Nr,GtNm1,'GtNm1 (THERMODYNAMICS)',myThid)
717           CALL DEBUG_STATS_RL(Nr,GsNm1,'GsNm1 (THERMODYNAMICS)',myThid)
718    #ifdef ALLOW_PTRACERS
719           IF ( usePTRACERS ) THEN
720             CALL PTRACERS_DEBUG(myThid)
721           ENDIF
722    #endif /* ALLOW_PTRACERS */
723          ENDIF
724    #endif
725    
726        RETURN        RETURN
727        END        END
Legend:
Removed from v.1.5  
changed lines
  Added in v.1.25
  ViewVC Help
Powered by ViewVC 1.1.22