/[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.4 by adcroft, Mon Sep 10 01:22:48 2001 UTC revision 1.102 by heimbach, Tue Apr 4 14:52:43 2006 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    
4    #include "PACKAGES_CONFIG.h"
5  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
6    
7    #ifdef ALLOW_AUTODIFF_TAMC
8    # ifdef ALLOW_GMREDI
9    #  include "GMREDI_OPTIONS.h"
10    # endif
11    # ifdef ALLOW_KPP
12    #  include "KPP_OPTIONS.h"
13    # endif
14    #endif /* ALLOW_AUTODIFF_TAMC */
15    
16    CBOP
17    C     !ROUTINE: THERMODYNAMICS
18    C     !INTERFACE:
19        SUBROUTINE THERMODYNAMICS(myTime, myIter, myThid)        SUBROUTINE THERMODYNAMICS(myTime, myIter, myThid)
20  C     /==========================================================\  C     !DESCRIPTION: \bv
21  C     | SUBROUTINE THERMODYNAMICS                                |  C     *==========================================================*
22  C     | o Controlling routine for the prognostic part of the     |  C     | SUBROUTINE THERMODYNAMICS
23  C     |   thermo-dynamics.                                       |  C     | o Controlling routine for the prognostic part of the
24  C     |==========================================================|  C     |   thermo-dynamics.
25  C     \==========================================================/  C     *===========================================================
26        IMPLICIT NONE  C     | The algorithm...
27    C     |
28    C     | "Correction Step"
29    C     | =================
30    C     | Here we update the horizontal velocities with the surface
31    C     | pressure such that the resulting flow is either consistent
32    C     | with the free-surface evolution or the rigid-lid:
33    C     |   U[n] = U* + dt x d/dx P
34    C     |   V[n] = V* + dt x d/dy P
35    C     |
36    C     | "Calculation of Gs"
37    C     | ===================
38    C     | This is where all the accelerations and tendencies (ie.
39    C     | physics, parameterizations etc...) are calculated
40    C     |   rho = rho ( theta[n], salt[n] )
41    C     |   b   = b(rho, theta)
42    C     |   K31 = K31 ( rho )
43    C     |   Gu[n] = Gu( u[n], v[n], wVel, b, ... )
44    C     |   Gv[n] = Gv( u[n], v[n], wVel, b, ... )
45    C     |   Gt[n] = Gt( theta[n], u[n], v[n], wVel, K31, ... )
46    C     |   Gs[n] = Gs( salt[n], u[n], v[n], wVel, K31, ... )
47    C     |
48    C     | "Time-stepping" or "Prediction"
49    C     | ================================
50    C     | The models variables are stepped forward with the appropriate
51    C     | time-stepping scheme (currently we use Adams-Bashforth II)
52    C     | - For momentum, the result is always *only* a "prediction"
53    C     | in that the flow may be divergent and will be "corrected"
54    C     | later with a surface pressure gradient.
55    C     | - Normally for tracers the result is the new field at time
56    C     | level [n+1} *BUT* in the case of implicit diffusion the result
57    C     | is also *only* a prediction.
58    C     | - We denote "predictors" with an asterisk (*).
59    C     |   U* = U[n] + dt x ( 3/2 Gu[n] - 1/2 Gu[n-1] )
60    C     |   V* = V[n] + dt x ( 3/2 Gv[n] - 1/2 Gv[n-1] )
61    C     |   theta[n+1] = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
62    C     |   salt[n+1] = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
63    C     | With implicit diffusion:
64    C     |   theta* = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
65    C     |   salt* = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
66    C     |   (1 + dt * K * d_zz) theta[n] = theta*
67    C     |   (1 + dt * K * d_zz) salt[n] = salt*
68    C     |
69    C     *==========================================================*
70    C     \ev
71    
72    C     !USES:
73          IMPLICIT NONE
74  C     == Global variables ===  C     == Global variables ===
75  #include "SIZE.h"  #include "SIZE.h"
76  #include "EEPARAMS.h"  #include "EEPARAMS.h"
77  #include "PARAMS.h"  #include "PARAMS.h"
78  #include "DYNVARS.h"  #include "DYNVARS.h"
79  #include "GRID.h"  #include "GRID.h"
80    #ifdef ALLOW_GENERIC_ADVDIFF
81  #include "GAD.h"  #include "GAD.h"
82  #ifdef ALLOW_PASSIVE_TRACER  #endif
83  #include "TR1.h"  #ifdef ALLOW_PTRACERS
84    #include "PTRACERS_SIZE.h"
85    #include "PTRACERS.h"
86    #endif
87    #ifdef ALLOW_TIMEAVE
88    #include "TIMEAVE_STATV.h"
89  #endif  #endif
90    
91  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
92  # include "tamc.h"  # include "tamc.h"
93  # include "tamc_keys.h"  # include "tamc_keys.h"
94  # include "FFIELDS.h"  # include "FFIELDS.h"
95    # include "EOS.h"
96  # ifdef ALLOW_KPP  # ifdef ALLOW_KPP
97  #  include "KPP.h"  #  include "KPP.h"
98  # endif  # endif
99  # ifdef ALLOW_GMREDI  # ifdef ALLOW_GMREDI
100  #  include "GMREDI.h"  #  include "GMREDI.h"
101  # endif  # endif
102    # ifdef ALLOW_EBM
103    #  include "EBM.h"
104    # endif
105  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
106    
107  #ifdef ALLOW_TIMEAVE  C     !INPUT/OUTPUT PARAMETERS:
 #include "TIMEAVE_STATV.h"  
 #endif  
   
108  C     == Routine arguments ==  C     == Routine arguments ==
109  C     myTime - Current time in simulation  C     myTime - Current time in simulation
110  C     myIter - Current iteration number in simulation  C     myIter - Current iteration number in simulation
# Line 47  C     myThid - Thread number for this in Line 113  C     myThid - Thread number for this in
113        INTEGER myIter        INTEGER myIter
114        INTEGER myThid        INTEGER myThid
115    
116    #ifdef ALLOW_GENERIC_ADVDIFF
117    C     !LOCAL VARIABLES:
118  C     == Local variables  C     == Local variables
119  C     xA, yA                 - Per block temporaries holding face areas  C     xA, yA                 - Per block temporaries holding face areas
120  C     uTrans, vTrans, rTrans - Per block temporaries holding flow  C     uTrans, vTrans, rTrans - Per block temporaries holding flow
# Line 54  C                              transport Line 122  C                              transport
122  C                              o uTrans: Zonal transport  C                              o uTrans: Zonal transport
123  C                              o vTrans: Meridional transport  C                              o vTrans: Meridional transport
124  C                              o rTrans: Vertical transport  C                              o rTrans: Vertical transport
125    C     rTransKp1                o vertical volume transp. at interface k+1
126  C     maskUp                   o maskUp: land/water mask for W points  C     maskUp                   o maskUp: land/water mask for W points
127  C     fVer[STUV]               o fVer: Vertical flux term - note fVer  C     fVer[STUV]               o fVer: Vertical flux term - note fVer
128  C                                      is "pipelined" in the vertical  C                                      is "pipelined" in the vertical
129  C                                      so we need an fVer for each  C                                      so we need an fVer for each
130  C                                      variable.  C                                      variable.
131  C     rhoK, rhoKM1   - Density at current level, and level above  C     kappaRT,       - Total diffusion in vertical at level k, for T and S
132  C     phiHyd         - Hydrostatic part of the potential phiHydi.  C     kappaRS          (background + spatially varying, isopycnal term).
133  C                      In z coords phiHydiHyd is the hydrostatic  C     kappaRTr       - Total diffusion in vertical at level k,
134  C                      Potential (=pressure/rho0) anomaly  C                      for each passive Tracer
135  C                      In p coords phiHydiHyd is the geopotential  C     kappaRk        - Total diffusion in vertical, all levels, 1 tracer  
136  C                      surface height anomaly.  C     useVariableK   = T when vertical diffusion is not constant
 C     phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean)  
 C     phiSurfY             or geopotentiel (atmos) in X and Y direction  
 C     KappaRT,       - Total diffusion in vertical for T and S.  
 C     KappaRS          (background + spatially varying, isopycnal term).  
137  C     iMin, iMax     - Ranges and sub-block indices on which calculations  C     iMin, iMax     - Ranges and sub-block indices on which calculations
138  C     jMin, jMax       are applied.  C     jMin, jMax       are applied.
139  C     bi, bj  C     bi, bj
# Line 80  C                      index into fVerTe Line 145  C                      index into fVerTe
145        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
146        _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
147        _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
148          _RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
149        _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
150        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
151        _RL fVerS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
152        _RL fVerTr1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL kappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
153        _RL phiHyd  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL kappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
154        _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)  #ifdef ALLOW_PTRACERS
155        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fVerP   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2,PTRACERS_num)
156        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL kappaRTr(1-Olx:sNx+Olx,1-Oly:sNy+Oly,PTRACERS_num)
157        _RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  #endif
158        _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL kappaRk (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
       _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)  
       _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
       _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
       _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
   
 C This is currently used by IVDC and Diagnostics  
       _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
   
159        INTEGER iMin, iMax        INTEGER iMin, iMax
160        INTEGER jMin, jMax        INTEGER jMin, jMax
161        INTEGER bi, bj        INTEGER bi, bj
162        INTEGER i, j        INTEGER i, j
163        INTEGER k, km1, kup, kDown        INTEGER k, km1, kup, kDown
164    #ifdef ALLOW_ADAMSBASHFORTH_3
165          INTEGER iterNb, m1, m2
166    #endif
167    #ifdef ALLOW_TIMEAVE
168          LOGICAL useVariableK
169    #endif
170    #ifdef ALLOW_PTRACERS
171          INTEGER iTracer, ip
172    #endif
173    
174  Cjmc : add for phiHyd output <- but not working if multi tile per CPU  CEOP
175  c     CHARACTER*(MAX_LEN_MBUF) suff  
176  c     LOGICAL  DIFFERENT_MULTIPLE  #ifdef ALLOW_DEBUG
177  c     EXTERNAL DIFFERENT_MULTIPLE           IF ( debugLevel .GE. debLevB )
178  Cjmc(end)       &    CALL DEBUG_ENTER('THERMODYNAMICS',myThid)
179    #endif
 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---  
180    
181  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
182  C--   dummy statement to end declaration part  C--   dummy statement to end declaration part
183        ikey = 1        ikey = 1
184          itdkey = 1
185  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
186    
 C--   Set up work arrays with valid (i.e. not NaN) values  
 C     These inital values do not alter the numerical results. They  
 C     just ensure that all memory references are to valid floating  
 C     point numbers. This prevents spurious hardware signals due to  
 C     uninitialised but inert locations.  
       DO j=1-OLy,sNy+OLy  
        DO i=1-OLx,sNx+OLx  
         xA(i,j)      = 0. _d 0  
         yA(i,j)      = 0. _d 0  
         uTrans(i,j)  = 0. _d 0  
         vTrans(i,j)  = 0. _d 0  
         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  
   
   
187  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
188  C--   HPF directive to help TAMC  C--   HPF directive to help TAMC
189  CHPF$ INDEPENDENT  CHPF$ INDEPENDENT
# Line 194  CHPF$ INDEPENDENT Line 194  CHPF$ INDEPENDENT
194  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
195  C--    HPF directive to help TAMC  C--    HPF directive to help TAMC
196  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS
197  CHPF$&                  ,phiHyd,utrans,vtrans,xA,yA  CHPF$&                  ,utrans,vtrans,xA,yA
198  CHPF$&                  ,KappaRT,KappaRS  CHPF$&                  ,kappaRT,kappaRS
199  CHPF$&                  )  CHPF$&                  )
200    # ifdef ALLOW_PTRACERS
201    CHPF$  INDEPENDENT, NEW (fVerP,kappaRTr)
202    # endif
203  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
204    
205         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
# Line 204  CHPF$&                  ) Line 207  CHPF$&                  )
207  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
208            act1 = bi - myBxLo(myThid)            act1 = bi - myBxLo(myThid)
209            max1 = myBxHi(myThid) - myBxLo(myThid) + 1            max1 = myBxHi(myThid) - myBxLo(myThid) + 1
   
210            act2 = bj - myByLo(myThid)            act2 = bj - myByLo(myThid)
211            max2 = myByHi(myThid) - myByLo(myThid) + 1            max2 = myByHi(myThid) - myByLo(myThid) + 1
   
212            act3 = myThid - 1            act3 = myThid - 1
213            max3 = nTx*nTy            max3 = nTx*nTy
   
214            act4 = ikey_dynamics - 1            act4 = ikey_dynamics - 1
215              itdkey = (act1 + 1) + act2*max1
           ikey = (act1 + 1) + act2*max1  
216       &                      + act3*max1*max2       &                      + act3*max1*max2
217       &                      + act4*max1*max2*max3       &                      + act4*max1*max2*max3
218  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
219    
220  C--     Set up work arrays that need valid initial values  C--   Set up work arrays with valid (i.e. not NaN) values
221    C     These inital values do not alter the numerical results. They
222    C     just ensure that all memory references are to valid floating
223    C     point numbers. This prevents spurious hardware signals due to
224    C     uninitialised but inert locations.
225    
226          DO j=1-OLy,sNy+OLy          DO j=1-OLy,sNy+OLy
227           DO i=1-OLx,sNx+OLx           DO i=1-OLx,sNx+OLx
228              xA(i,j)        = 0. _d 0
229              yA(i,j)        = 0. _d 0
230              uTrans(i,j)    = 0. _d 0
231              vTrans(i,j)    = 0. _d 0
232            rTrans (i,j)   = 0. _d 0            rTrans (i,j)   = 0. _d 0
233              rTransKp1(i,j) = 0. _d 0
234            fVerT  (i,j,1) = 0. _d 0            fVerT  (i,j,1) = 0. _d 0
235            fVerT  (i,j,2) = 0. _d 0            fVerT  (i,j,2) = 0. _d 0
236            fVerS  (i,j,1) = 0. _d 0            fVerS  (i,j,1) = 0. _d 0
237            fVerS  (i,j,2) = 0. _d 0            fVerS  (i,j,2) = 0. _d 0
238            fVerTr1(i,j,1) = 0. _d 0            kappaRT(i,j)   = 0. _d 0
239            fVerTr1(i,j,2) = 0. _d 0            kappaRS(i,j)   = 0. _d 0
240           ENDDO           ENDDO
241          ENDDO          ENDDO
242    
# Line 235  C--     Set up work arrays that need val Line 244  C--     Set up work arrays that need val
244           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
245            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
246  C This is currently also used by IVDC and Diagnostics  C This is currently also used by IVDC and Diagnostics
247             ConvectCount(i,j,k) = 0.             kappaRk(i,j,k)    = 0. _d 0
248             KappaRT(i,j,k) = 0. _d 0  C- tracer tendency needs to be set to zero (moved here from gad_calc_rhs):
249             KappaRS(i,j,k) = 0. _d 0             gT(i,j,k,bi,bj)   = 0. _d 0
250               gS(i,j,k,bi,bj)   = 0. _d 0
251            ENDDO            ENDDO
252           ENDDO           ENDDO
253          ENDDO          ENDDO
254    
255          iMin = 1-OLx+1  #ifdef ALLOW_PTRACERS
256          iMax = sNx+OLx          IF ( usePTRACERS ) THEN
257          jMin = 1-OLy+1           DO ip=1,PTRACERS_num
258          jMax = sNy+OLy             DO j=1-OLy,sNy+OLy
259                DO i=1-OLx,sNx+OLx
260                 fVerP  (i,j,1,ip) = 0. _d 0
261  #ifdef ALLOW_AUTODIFF_TAMC               fVerP  (i,j,2,ip) = 0. _d 0
262  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte               kappaRTr(i,j,ip)  = 0. _d 0
263  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte              ENDDO
264  #endif /* ALLOW_AUTODIFF_TAMC */             ENDDO
265             ENDDO
266  C--     Start of diagnostic loop  C-      set tracer tendency to zero:
267          DO k=Nr,1,-1           DO iTracer=1,PTRACERS_num
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 C? Patrick, is this formula correct now that we change the loop range?  
 C? Do we still need this?  
 cph kkey formula corrected.  
 cph Needed for rhok, rhokm1, in the case useGMREDI.  
          kkey = (ikey-1)*Nr + k  
 CADJ STORE rhokm1(:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte  
 CADJ STORE rhok  (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 C--       Integrate continuity vertically for vertical velocity  
           CALL INTEGRATE_FOR_W(  
      I                         bi, bj, k, uVel, vVel,  
      O                         wVel,  
      I                         myThid )  
   
 #ifdef    ALLOW_OBCS  
 #ifdef    ALLOW_NONHYDROSTATIC  
 C--       Apply OBC to W if in N-H mode  
           IF (useOBCS.AND.nonHydrostatic) THEN  
             CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )  
           ENDIF  
 #endif    /* ALLOW_NONHYDROSTATIC */  
 #endif    /* ALLOW_OBCS */  
   
 C--       Calculate gradients of potential density for isoneutral  
 C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)  
 c         IF ( k.GT.1 .AND. (useGMRedi.OR.ivdc_kappa.NE.0.) ) THEN  
           IF ( useGMRedi .OR. (k.GT.1 .AND. ivdc_kappa.NE.0.) ) THEN  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  
 CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
             CALL FIND_RHO(  
      I        bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,  
      I        theta, salt,  
      O        rhoK,  
      I        myThid )  
             IF (k.GT.1) THEN  
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE theta(:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  
 CADJ STORE salt (:,:,k-1,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
              CALL FIND_RHO(  
      I        bi, bj, iMin, iMax, jMin, jMax, k-1, k, eosType,  
      I        theta, salt,  
      O        rhoKm1,  
      I        myThid )  
             ENDIF  
             CALL GRAD_SIGMA(  
      I             bi, bj, iMin, iMax, jMin, jMax, k,  
      I             rhoK, rhoKm1, rhoK,  
      O             sigmaX, sigmaY, sigmaR,  
      I             myThid )  
           ENDIF  
   
 C--       Implicit Vertical Diffusion for Convection  
 c ==> should use sigmaR !!!  
           IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN  
             CALL CALC_IVDC(  
      I        bi, bj, iMin, iMax, jMin, jMax, k,  
      I        rhoKm1, rhoK,  
      U        ConvectCount, KappaRT, KappaRS,  
      I        myTime, myIter, myThid)  
           ENDIF  
   
 C--     end of diagnostic k loop (Nr:1)  
         ENDDO  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 cph avoids recomputation of integrate_for_w  
 CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 #ifdef  ALLOW_OBCS  
 C--     Calculate future values on open boundaries  
         IF (useOBCS) THEN  
           CALL OBCS_CALC( bi, bj, myTime+deltaT,  
      I            uVel, vVel, wVel, theta, salt,  
      I            myThid )  
         ENDIF  
 #endif  /* ALLOW_OBCS */  
   
 C--     Determines forcing terms based on external fields  
 C       relaxation terms, etc.  
         CALL EXTERNAL_FORCING_SURF(  
      I             bi, bj, iMin, iMax, jMin, jMax,  
      I             myThid )  
 #ifdef ALLOW_AUTODIFF_TAMC  
 cph needed for KPP  
 CADJ STORE surfacetendencyU(:,:,bi,bj)  
 CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE surfacetendencyV(:,:,bi,bj)  
 CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE surfacetendencyS(:,:,bi,bj)  
 CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  
 CADJ STORE surfacetendencyT(:,:,bi,bj)  
 CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
   
 #ifdef  ALLOW_GMREDI  
   
 #ifdef ALLOW_AUTODIFF_TAMC  
 CADJ STORE sigmaX(:,:,:) = comlev1, key=ikey, byte=isbyte  
 CADJ STORE sigmaY(:,:,:) = comlev1, key=ikey, byte=isbyte  
 CADJ STORE sigmaR(:,:,:) = comlev1, key=ikey, byte=isbyte  
 #endif /* ALLOW_AUTODIFF_TAMC */  
 C--     Calculate iso-neutral slopes for the GM/Redi parameterisation  
         IF (useGMRedi) THEN  
268            DO k=1,Nr            DO k=1,Nr
269              CALL GMREDI_CALC_TENSOR(             DO j=1-OLy,sNy+OLy
270       I             bi, bj, iMin, iMax, jMin, jMax, k,              DO i=1-OLx,sNx+OLx
271       I             sigmaX, sigmaY, sigmaR,               gPTr(i,j,k,bi,bj,itracer) = 0. _d 0
272       I             myThid )              ENDDO
273            ENDDO             ENDDO
 #ifdef ALLOW_AUTODIFF_TAMC  
         ELSE  
           DO k=1, Nr  
             CALL GMREDI_CALC_TENSOR_DUMMY(  
      I             bi, bj, iMin, iMax, jMin, jMax, k,  
      I             sigmaX, sigmaY, sigmaR,  
      I             myThid )  
274            ENDDO            ENDDO
275  #endif /* ALLOW_AUTODIFF_TAMC */           ENDDO
276          ENDIF          ENDIF
277    #endif
278    
279  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_ADAMSBASHFORTH_3
280  CADJ STORE Kwx(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte  C-      Apply AB on T,S :
281  CADJ STORE Kwy(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte          iterNb = myIter
282  CADJ STORE Kwz(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte          IF (staggerTimeStep) iterNb = myIter - 1
283  #endif /* ALLOW_AUTODIFF_TAMC */          m1 = 1 + MOD(iterNb+1,2)
284            m2 = 1 + MOD( iterNb ,2)
285  #endif  /* ALLOW_GMREDI */  C       compute T^n+1/2 (stored in gtNm) extrapolating T forward in time
286            IF ( AdamsBashforth_T ) CALL ADAMS_BASHFORTH3(
287  #ifdef  ALLOW_KPP       I                                  bi, bj, 0,
288  C--     Compute KPP mixing coefficients       U                                  theta, gtNm,
289          IF (useKPP) THEN       I                                  tempStartAB, iterNb, myThid )
290            CALL KPP_CALC(  C       compute S^n+1/2 (stored in gsNm) extrapolating S forward in time
291       I                  bi, bj, myTime, myThid )          IF ( AdamsBashforth_S ) CALL ADAMS_BASHFORTH3(
292  #ifdef ALLOW_AUTODIFF_TAMC       I                                  bi, bj, 0,
293          ELSE       U                                  salt, gsNm,
294            CALL KPP_CALC_DUMMY(       I                                  saltStartAB, iterNb, myThid )
295       I                  bi, bj, myTime, myThid )  #endif /* ALLOW_ADAMSBASHFORTH_3 */
296  #endif /* ALLOW_AUTODIFF_TAMC */  
297          ENDIF  c       iMin = 1-OLx
298    c       iMax = sNx+OLx
299    c       jMin = 1-OLy
300    c       jMax = sNy+OLy
301    
302  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
303  CADJ STORE KPPghat   (:,:,:,bi,bj)  cph avoids recomputation of integrate_for_w
304  CADJ &   , KPPviscAz (:,:,:,bi,bj)  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
 CADJ &   , KPPdiffKzT(:,:,:,bi,bj)  
 CADJ &   , KPPdiffKzS(:,:,:,bi,bj)  
 CADJ &   , KPPfrac   (:,:  ,bi,bj)  
 CADJ &                 = comlev1_bibj, key=ikey, byte=isbyte  
305  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
306    
307  #endif  /* ALLOW_KPP */  C--     Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h
308    C--     MOST of THERMODYNAMICS will be disabled
309    #ifndef SINGLE_LAYER_MODE
310    
311  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
312  CADJ STORE KappaRT(:,:,:)     = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
313  CADJ STORE KappaRS(:,:,:)     = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
314  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
315  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
316  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  cph-test
317  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  cphCADJ STORE gtnm1(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
318  #ifdef ALLOW_PASSIVE_TRACER  cphCADJ STORE gsnm1(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
 CADJ STORE tr1  (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
 #endif  
319  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
320    
321  #ifdef ALLOW_AIM  #ifndef DISABLE_MULTIDIM_ADVECTION
 C       AIM - atmospheric intermediate model, physics package code.  
 C note(jmc) : phiHyd=0 at this point but is not really used in Molteni Physics  
         IF ( useAIM ) THEN  
          CALL TIMER_START('AIM_DO_ATMOS_PHYS      [DYNAMICS]', myThid)  
          CALL AIM_DO_ATMOS_PHYSICS( phiHyd, bi, bj, myTime, myThid )  
          CALL TIMER_STOP ('AIM_DO_ATMOS_PHYS      [DYNAMICS]', myThid)  
         ENDIF  
 #endif /* ALLOW_AIM */  
   
322  C--     Some advection schemes are better calculated using a multi-dimensional  C--     Some advection schemes are better calculated using a multi-dimensional
323  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.
324          IF (multiDimAdvection) THEN  C
325           IF (tempStepping .AND.  C The CPP flag DISABLE_MULTIDIM_ADVECTION is currently unset in GAD_OPTIONS.h
326       &       tempAdvScheme.NE.ENUM_CENTERED_2ND .AND.  C The default is to use multi-dimensinal advection for non-linear advection
327       &       tempAdvScheme.NE.ENUM_UPWIND_3RD .AND.  C schemes. However, for the sake of efficiency of the adjoint it is necessary
328       &       tempAdvScheme.NE.ENUM_CENTERED_4TH )  C to be able to exclude this scheme to avoid excessive storage and
329       &   CALL GAD_ADVECTION(bi,bj,tempAdvScheme,GAD_TEMPERATURE,theta,  C recomputation. It *is* differentiable, if you need it.
330       U                      gT,  C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to
331       I                      myTime,myIter,myThid)  C disable this section of code.
332           IF (saltStepping .AND.          IF (tempMultiDimAdvec) THEN
333       &       saltAdvScheme.NE.ENUM_CENTERED_2ND .AND.  #ifdef ALLOW_DEBUG
334       &       saltAdvScheme.NE.ENUM_UPWIND_3RD .AND.            IF ( debugLevel .GE. debLevB )
335       &       saltAdvScheme.NE.ENUM_CENTERED_4TH )       &     CALL DEBUG_CALL('GAD_ADVECTION',myThid)
336       &   CALL GAD_ADVECTION(bi,bj,saltAdvScheme,GAD_SALINITY,salt,  #endif
337       U                      gS,            CALL GAD_ADVECTION(
338       I                      myTime,myIter,myThid)       I             tempImplVertAdv, tempAdvScheme, tempVertAdvScheme,
339         I             GAD_TEMPERATURE,
340         I             uVel, vVel, wVel, theta,
341         O             gT,
342         I             bi,bj,myTime,myIter,myThid)
343            ENDIF
344            IF (saltMultiDimAdvec) THEN
345    #ifdef ALLOW_DEBUG
346              IF ( debugLevel .GE. debLevB )
347         &     CALL DEBUG_CALL('GAD_ADVECTION',myThid)
348    #endif
349              CALL GAD_ADVECTION(
350         I             saltImplVertAdv, saltAdvScheme, saltVertAdvScheme,
351         I             GAD_SALINITY,
352         I             uVel, vVel, wVel, salt,
353         O             gS,
354         I             bi,bj,myTime,myIter,myThid)
355          ENDIF          ENDIF
356    
357    C Since passive tracers are configurable separately from T,S we
358    C call the multi-dimensional method for PTRACERS regardless
359    C of whether multiDimAdvection is set or not.
360    #ifdef ALLOW_PTRACERS
361            IF ( usePTRACERS ) THEN
362    #ifdef ALLOW_DEBUG
363              IF ( debugLevel .GE. debLevB )
364         &     CALL DEBUG_CALL('PTRACERS_ADVECTION',myThid)
365    #endif
366             CALL PTRACERS_ADVECTION( bi,bj,myIter,myTime,myThid )
367            ENDIF
368    #endif /* ALLOW_PTRACERS */
369    #endif /* DISABLE_MULTIDIM_ADVECTION */
370    
371    #ifdef ALLOW_DEBUG
372            IF ( debugLevel .GE. debLevB )
373         &    CALL DEBUG_MSG('ENTERING DOWNWARD K LOOP',myThid)
374    #endif
375    
376  C--     Start of thermodynamics loop  C--     Start of thermodynamics loop
377          DO k=Nr,1,-1          DO k=Nr,1,-1
378  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
379  C? Patrick Is this formula correct?  C? Patrick Is this formula correct?
380  cph Yes, but I rewrote it.  cph Yes, but I rewrote it.
381  cph Also, the KappaR? need the index and subscript k!  cph Also, the kappaR? need the index and subscript k!
382           kkey = (ikey-1)*Nr + k           kkey = (itdkey-1)*Nr + k
383  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
384    
385  C--       km1    Points to level above k (=k-1)  C--       km1    Points to level above k (=k-1)
# Line 477  C--       kDown  Cycles through 2,1 to p Line 395  C--       kDown  Cycles through 2,1 to p
395            jMin = 1-OLy            jMin = 1-OLy
396            jMax = sNy+OLy            jMax = sNy+OLy
397    
398              IF (k.EQ.Nr) THEN
399               DO j=1-Oly,sNy+Oly
400                DO i=1-Olx,sNx+Olx
401                 rTransKp1(i,j) = 0. _d 0
402                ENDDO
403               ENDDO
404              ELSE
405               DO j=1-Oly,sNy+Oly
406                DO i=1-Olx,sNx+Olx
407                 rTransKp1(i,j) = rTrans(i,j)
408                ENDDO
409               ENDDO
410              ENDIF
411    #ifdef ALLOW_AUTODIFF_TAMC
412    CADJ STORE rTransKp1(:,:) = comlev1_bibj_k, key=kkey, byte=isbyte
413    #endif
414    
415  C--       Get temporary terms used by tendency routines  C--       Get temporary terms used by tendency routines
416            CALL CALC_COMMON_FACTORS (            CALL CALC_COMMON_FACTORS (
417       I         bi,bj,iMin,iMax,jMin,jMax,k,       I         bi,bj,iMin,iMax,jMin,jMax,k,
418       O         xA,yA,uTrans,vTrans,rTrans,maskUp,       O         xA,yA,uTrans,vTrans,rTrans,maskUp,
419       I         myThid)       I         myThid)
420    
421  #ifdef ALLOW_AUTODIFF_TAMC            IF (k.EQ.1) THEN
422  CADJ STORE KappaRT(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte  C- Surface interface :
423  CADJ STORE KappaRS(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte             DO j=1-Oly,sNy+Oly
424                DO i=1-Olx,sNx+Olx
425                 rTrans(i,j) = 0.
426                ENDDO
427               ENDDO
428              ELSE
429    C- Interior interface :
430               DO j=1-Oly,sNy+Oly
431                DO i=1-Olx,sNx+Olx
432                 rTrans(i,j) = rTrans(i,j)*maskC(i,j,k-1,bi,bj)
433                ENDDO
434               ENDDO
435              ENDIF
436    
437    #ifdef ALLOW_GMREDI
438    
439    C--   Residual transp = Bolus transp + Eulerian transp
440              IF (useGMRedi) THEN
441                CALL GMREDI_CALC_UVFLOW(
442         &            uTrans, vTrans, bi, bj, k, myThid)
443                IF (K.GE.2) CALL GMREDI_CALC_WFLOW(
444         &                    rTrans, bi, bj, k, myThid)
445              ENDIF
446    
447    #ifdef ALLOW_AUTODIFF_TAMC
448    CADJ STORE rTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
449    #ifdef GM_BOLUS_ADVEC
450    CADJ STORE uTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
451    CADJ STORE vTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
452    #endif
453  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
454    
455    #endif /* ALLOW_GMREDI */
456    
457  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL  #ifdef  INCLUDE_CALC_DIFFUSIVITY_CALL
458  C--      Calculate the total vertical diffusivity  C--      Calculate the total vertical diffusivity
459           CALL CALC_DIFFUSIVITY(            IF ( .NOT.implicitDiffusion ) THEN
460       I        bi,bj,iMin,iMax,jMin,jMax,k,              CALL CALC_DIFFUSIVITY(
461       I        maskUp,       I          bi,bj,iMin,iMax,jMin,jMax,k,
462       O        KappaRT,KappaRS,       I          maskUp,
463       I        myThid)       O          kappaRT,kappaRS,
464         I          myThid)
465              ENDIF
466    # ifdef ALLOW_AUTODIFF_TAMC
467    CADJ STORE kappaRT(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
468    CADJ STORE kappaRS(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
469    # endif /* ALLOW_AUTODIFF_TAMC */
470  #endif  #endif
471    
472            iMin = 1-OLx+2            iMin = 1-OLx+2
# Line 503  C--      Calculate the total vertical di Line 475  C--      Calculate the total vertical di
475            jMax = sNy+OLy-1            jMax = sNy+OLy-1
476    
477  C--      Calculate active tracer tendencies (gT,gS,...)  C--      Calculate active tracer tendencies (gT,gS,...)
478  C        and step forward storing result in gTnm1, gSnm1, etc.  C        and step forward storing result in gT, gS, etc.
479    C--
480    # ifdef ALLOW_AUTODIFF_TAMC
481    #  if (defined (NONLIN_FRSURF) && defined (ALLOW_GMREDI))
482    CADJ STORE kux(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
483    CADJ STORE kuz(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
484    CADJ STORE kvy(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
485    CADJ STORE kvz(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
486    #  endif
487    # endif /* ALLOW_AUTODIFF_TAMC */
488           IF ( tempStepping ) THEN           IF ( tempStepping ) THEN
489    #ifdef ALLOW_AUTODIFF_TAMC
490    CADJ STORE gTnm1(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
491    # ifdef NONLIN_FRSURF
492    CADJ STORE gT(:,:,k,bi,bj)    = comlev1_bibj_k, key=kkey, byte=isbyte
493    # endif
494    #endif /* ALLOW_AUTODIFF_TAMC */
495             CALL CALC_GT(             CALL CALC_GT(
496       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
497       I         xA,yA,uTrans,vTrans,rTrans,maskUp,       I         xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,
498       I         KappaRT,       I         kappaRT,
499       U         fVerT,       U         fVerT,
500       I         myTime, myThid)       I         myTime,myIter,myThid)
501    #ifdef ALLOW_ADAMSBASHFORTH_3
502              IF ( AdamsBashforth_T ) THEN
503               CALL TIMESTEP_TRACER(
504         I         bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,
505         I         gtNm(1-Olx,1-Oly,1,1,1,m2),
506         U         gT,
507         I         myIter, myThid)
508              ELSE
509    #endif
510             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
511       I         bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,       I         bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,
512       I         theta, gT,       I         theta,
513       U         gTnm1,       U         gT,
514       I         myIter, myThid)       I         myIter, myThid)
515    #ifdef ALLOW_ADAMSBASHFORTH_3
516              ENDIF
517    #endif
518           ENDIF           ENDIF
519    
520           IF ( saltStepping ) THEN           IF ( saltStepping ) THEN
521    #ifdef ALLOW_AUTODIFF_TAMC
522    CADJ STORE gSnm1(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte
523    # ifdef NONLIN_FRSURF
524    CADJ STORE gS(:,:,k,bi,bj)    = comlev1_bibj_k, key=kkey, byte=isbyte
525    # endif
526    #endif /* ALLOW_AUTODIFF_TAMC */
527             CALL CALC_GS(             CALL CALC_GS(
528       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
529       I         xA,yA,uTrans,vTrans,rTrans,maskUp,       I         xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,
530       I         KappaRS,       I         kappaRS,
531       U         fVerS,       U         fVerS,
532       I         myTime, myThid)       I         myTime,myIter,myThid)
533    #ifdef ALLOW_ADAMSBASHFORTH_3
534              IF ( AdamsBashforth_S ) THEN
535             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
536       I         bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,       I         bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,
537       I         salt, gS,       I         gsNm(1-Olx,1-Oly,1,1,1,m2),
538       U         gSnm1,       U         gS,
539       I         myIter, myThid)       I         myIter, myThid)
540           ENDIF            ELSE
541  #ifdef ALLOW_PASSIVE_TRACER  #endif
          IF ( tr1Stepping ) THEN  
            CALL CALC_GTR1(  
      I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,  
      I         xA,yA,uTrans,vTrans,rTrans,maskUp,  
      I         KappaRT,  
      U         fVerTr1,  
      I         myTime, myThid)  
542             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
543       I         bi,bj,iMin,iMax,jMin,jMax,k,tracerAdvScheme,       I         bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,
544       I         Tr1, gTr1,       I         salt,
545       U         gTr1NM1,       U         gS,
546       I         myIter, myThid)       I         myIter, myThid)
547           ENDIF  #ifdef ALLOW_ADAMSBASHFORTH_3
548              ENDIF
549  #endif  #endif
550             ENDIF
551    
552    #ifdef ALLOW_PTRACERS
553             IF ( usePTRACERS ) THEN
554               IF ( .NOT.implicitDiffusion ) THEN
555                 CALL PTRACERS_CALC_DIFF(
556         I            bi,bj,iMin,iMax,jMin,jMax,k,
557         I            maskUp,
558         O            kappaRTr,
559         I            myThid)
560               ENDIF
561    # ifdef ALLOW_AUTODIFF_TAMC
562    CADJ STORE kappaRTr(:,:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
563    # endif /* ALLOW_AUTODIFF_TAMC */
564               CALL PTRACERS_INTEGRATE(
565         I         bi,bj,k,
566         I         xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,
567         U         fVerP,
568         I         kappaRTr,
569         I         myIter,myTime,myThid)
570             ENDIF
571    #endif /* ALLOW_PTRACERS */
572    
573  #ifdef   ALLOW_OBCS  #ifdef   ALLOW_OBCS
574  C--      Apply open boundary conditions  C--      Apply open boundary conditions
575           IF (useOBCS) THEN           IF (useOBCS) THEN
576             CALL OBCS_APPLY_TS( bi, bj, k, gTnm1, gSnm1, myThid )             CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )
577           END IF           END IF
578  #endif   /* ALLOW_OBCS */  #endif   /* ALLOW_OBCS */
579    
580  C--      Freeze water  C--      Freeze water
581           IF (allowFreezing) THEN  C  this bit of code is left here for backward compatibility.
582    C  freezing at surface level has been moved to FORWARD_STEP
583             IF ( useOldFreezing .AND. .NOT. useSEAICE
584         &       .AND. .NOT.(useThSIce.AND.k.EQ.1) ) THEN
585  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
586  CADJ STORE gTNm1(:,:,k,bi,bj) = comlev1_bibj_k  CADJ STORE gT(:,:,k,bi,bj) = comlev1_bibj_k
587  CADJ &   , key = kkey, byte = isbyte  CADJ &   , key = kkey, byte = isbyte
588  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
589              CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )              CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )
590           END IF           ENDIF
591    
592  C--     end of thermodynamic k loop (Nr:1)  C--     end of thermodynamic k loop (Nr:1)
593          ENDDO          ENDDO
594    
595    C       All explicit advection/diffusion/sources should now be
596  #ifdef ALLOW_AUTODIFF_TAMC  C       done. The updated tracer field is in gPtr. Accumalate
597  C? Patrick? What about this one?  C       explicit tendency and also reset gPtr to initial tracer
598  cph Keys iikey and idkey don't seem to be needed  C       field for implicit matrix calculation
599  cph since storing occurs on different tape for each  
600  cph impldiff call anyways.  #ifdef ALLOW_MATRIX
601  cph Thus, common block comlev1_impl isn't needed either.          IF (useMATRIX)
602  cph Storing below needed in the case useGMREDI.       &    CALL MATRIX_STORE_TENDENCY_EXP(bi,bj, myTime,myIter,myThid)
603          iikey = (ikey-1)*maximpl  #endif
604  #endif /* ALLOW_AUTODIFF_TAMC */  
605            iMin = 1
606  C--     Implicit diffusion          iMax = sNx
607          IF (implicitDiffusion) THEN          jMin = 1
608            jMax = sNy
609           IF (tempStepping) THEN  
610    C--     Implicit vertical advection & diffusion
611            IF ( tempStepping .AND. implicitDiffusion ) THEN
612              CALL CALC_3D_DIFFUSIVITY(
613         I         bi,bj,iMin,iMax,jMin,jMax,
614         I         GAD_TEMPERATURE, useGMredi, useKPP,
615         O         kappaRk,
616         I         myThid)
617            ENDIF
618    #ifdef INCLUDE_IMPLVERTADV_CODE
619            IF ( tempImplVertAdv ) THEN
620              CALL GAD_IMPLICIT_R(
621         I         tempImplVertAdv, tempAdvScheme, GAD_TEMPERATURE,
622         I         kappaRk, wVel, theta,
623         U         gT,
624         I         bi, bj, myTime, myIter, myThid )
625            ELSEIF ( tempStepping .AND. implicitDiffusion ) THEN
626    #else /* INCLUDE_IMPLVERTADV_CODE */
627            IF     ( tempStepping .AND. implicitDiffusion ) THEN
628    #endif /* INCLUDE_IMPLVERTADV_CODE */
629  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
630              idkey = iikey + 1  CADJ STORE kappaRk(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
631  CADJ STORE gTNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
632  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
633              CALL IMPLDIFF(            CALL IMPLDIFF(
634       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
635       I         deltaTtracer, KappaRT, recip_HFacC,       I         GAD_TEMPERATURE, kappaRk, recip_hFacC,
636       U         gTNm1,       U         gT,
637       I         myThid )       I         myThid )
638            ENDIF
639    
640    #ifdef ALLOW_TIMEAVE
641            useVariableK = useKPP .OR. usePP81 .OR. useMY82 .OR. useGGL90
642         &       .OR. useGMredi .OR. ivdc_kappa.NE.0.
643            IF (taveFreq.GT.0. .AND. useVariableK ) THEN
644             IF (implicitDiffusion) THEN
645               CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, gT, kappaRk,
646         I                        Nr, 3, deltaTclock, bi, bj, myThid)
647    c        ELSE
648    c          CALL TIMEAVE_CUMUL_DIF_1T(TdiffRtave, theta, kappaRT,
649    c    I                        Nr, 3, deltaTclock, bi, bj, myThid)
650           ENDIF           ENDIF
651            ENDIF
652    #endif /* ALLOW_TIMEAVE */
653    
654            IF ( saltStepping .AND. implicitDiffusion ) THEN
655              CALL CALC_3D_DIFFUSIVITY(
656         I         bi,bj,iMin,iMax,jMin,jMax,
657         I         GAD_SALINITY, useGMredi, useKPP,
658         O         kappaRk,
659         I         myThid)
660            ENDIF
661    
662           IF (saltStepping) THEN  #ifdef INCLUDE_IMPLVERTADV_CODE
663            IF ( saltImplVertAdv ) THEN
664              CALL GAD_IMPLICIT_R(
665         I         saltImplVertAdv, saltAdvScheme, GAD_SALINITY,
666         I         kappaRk, wVel, salt,
667         U         gS,
668         I         bi, bj, myTime, myIter, myThid )
669            ELSEIF ( saltStepping .AND. implicitDiffusion ) THEN
670    #else /* INCLUDE_IMPLVERTADV_CODE */
671            IF     ( saltStepping .AND. implicitDiffusion ) THEN
672    #endif /* INCLUDE_IMPLVERTADV_CODE */
673  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
674           idkey = iikey + 2  CADJ STORE kappaRk(:,:,:) = comlev1_bibj , key=itdkey, byte=isbyte
675  CADJ STORE gSNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
676  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
677              CALL IMPLDIFF(            CALL IMPLDIFF(
678       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
679       I         deltaTtracer, KappaRS, recip_HFacC,       I         GAD_SALINITY, kappaRk, recip_hFacC,
680       U         gSNm1,       U         gS,
681       I         myThid )       I         myThid )
682           ENDIF          ENDIF
683    
684  #ifdef ALLOW_PASSIVE_TRACER  #ifdef ALLOW_PTRACERS
685           IF (tr1Stepping) THEN          IF     ( usePTRACERS ) THEN
686  #ifdef ALLOW_AUTODIFF_TAMC  C--     Vertical advection/diffusion (implicit) for passive tracers
687  CADJ STORE gTr1Nm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte             CALL PTRACERS_IMPLICIT(
688  #endif /* ALLOW_AUTODIFF_TAMC */       U                             kappaRk,
689            CALL IMPLDIFF(       I                             bi, bj, myTime, myIter, myThid )
690       I      bi, bj, iMin, iMax, jMin, jMax,          ENDIF
691       I      deltaTtracer, KappaRT, recip_HFacC,  #endif /* ALLOW_PTRACERS */
      U      gTr1Nm1,  
      I      myThid )  
          ENDIF  
 #endif  
692    
693  #ifdef   ALLOW_OBCS  #ifdef   ALLOW_OBCS
694  C--      Apply open boundary conditions  C--      Apply open boundary conditions
695           IF (useOBCS) THEN          IF ( ( implicitDiffusion
696         &    .OR. tempImplVertAdv
697         &    .OR. saltImplVertAdv
698         &       ) .AND. useOBCS     ) THEN
699             DO K=1,Nr             DO K=1,Nr
700               CALL OBCS_APPLY_TS( bi, bj, k, gTnm1, gSnm1, myThid )               CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )
701             ENDDO             ENDDO
702           END IF          ENDIF
703  #endif   /* ALLOW_OBCS */  #endif   /* ALLOW_OBCS */
704    
705  C--     End If implicitDiffusion  #endif /* SINGLE_LAYER_MODE */
         ENDIF  
706    
707  Ccs-  C--   end bi,bj loops.
708         ENDDO         ENDDO
709        ENDDO        ENDDO
710    
711  #ifdef ALLOW_AIM  #ifdef ALLOW_DEBUG
712        IF ( useAIM ) THEN        If (debugMode) THEN
713         CALL AIM_AIM2DYN_EXCHANGES( myTime, myThid )         CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (THERMODYNAMICS)',myThid)
714        ENDIF         CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (THERMODYNAMICS)',myThid)
715         _EXCH_XYZ_R8(gTnm1,myThid)         CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (THERMODYNAMICS)',myThid)
716         _EXCH_XYZ_R8(gSnm1,myThid)         CALL DEBUG_STATS_RL(Nr,theta,'Theta (THERMODYNAMICS)',myThid)
717  #else         CALL DEBUG_STATS_RL(Nr,salt,'Salt (THERMODYNAMICS)',myThid)
718        IF (staggerTimeStep.AND.useCubedSphereExchange) THEN         CALL DEBUG_STATS_RL(Nr,gT,'Gt (THERMODYNAMICS)',myThid)
719         _EXCH_XYZ_R8(gTnm1,myThid)         CALL DEBUG_STATS_RL(Nr,gS,'Gs (THERMODYNAMICS)',myThid)
720         _EXCH_XYZ_R8(gSnm1,myThid)  #ifndef ALLOW_ADAMSBASHFORTH_3
721           CALL DEBUG_STATS_RL(Nr,gtNm1,'GtNm1 (THERMODYNAMICS)',myThid)
722           CALL DEBUG_STATS_RL(Nr,gsNm1,'GsNm1 (THERMODYNAMICS)',myThid)
723    #endif
724    #ifdef ALLOW_PTRACERS
725           IF ( usePTRACERS ) THEN
726             CALL PTRACERS_DEBUG(myThid)
727           ENDIF
728    #endif /* ALLOW_PTRACERS */
729        ENDIF        ENDIF
730  #endif /* ALLOW_AIM */  #endif /* ALLOW_DEBUG */
731    
732    #ifdef ALLOW_DEBUG
733             IF ( debugLevel .GE. debLevB )
734         &    CALL DEBUG_LEAVE('THERMODYNAMICS',myThid)
735    #endif
736    
737    #endif /* ALLOW_GENERIC_ADVDIFF */
738    
739        RETURN        RETURN
740        END        END

Legend:
Removed from v.1.4  
changed lines
  Added in v.1.102

  ViewVC Help
Powered by ViewVC 1.1.22