/[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.1 by adcroft, Fri Aug 3 19:06:11 2001 UTC revision 1.34 by heimbach, Fri Jan 10 19:06:05 2003 UTC
# Line 2  C $Header$ Line 2  C $Header$
2  C $Name$  C $Name$
3    
4  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
5    #ifdef ALLOW_AUTODIFF_TAMC
6    # ifdef ALLOW_GMREDI
7    #  include "GMREDI_OPTIONS.h"
8    # endif
9    # ifdef ALLOW_KPP
10    #  include "KPP_OPTIONS.h"
11    # endif
12    cswdice --- add ----
13    #ifdef ALLOW_THERM_SEAICE
14    #include "ICE.h"
15    #endif
16    cswdice ------
17    #endif /* ALLOW_AUTODIFF_TAMC */
18    
19    CBOP
20    C     !ROUTINE: THERMODYNAMICS
21    C     !INTERFACE:
22        SUBROUTINE THERMODYNAMICS(myTime, myIter, myThid)        SUBROUTINE THERMODYNAMICS(myTime, myIter, myThid)
23  C     /==========================================================\  C     !DESCRIPTION: \bv
24  C     | SUBROUTINE THERMODYNAMICS                                |  C     *==========================================================*
25  C     | o Controlling routine for the prognostic part of the     |  C     | SUBROUTINE THERMODYNAMICS                                
26  C     |   thermo-dynamics.                                       |  C     | o Controlling routine for the prognostic part of the      
27  C     |==========================================================|  C     |   thermo-dynamics.                                        
28  C     \==========================================================/  C     *===========================================================
29        IMPLICIT NONE  C     | The algorithm...
30    C     |
31    C     | "Correction Step"
32    C     | =================
33    C     | Here we update the horizontal velocities with the surface
34    C     | pressure such that the resulting flow is either consistent
35    C     | with the free-surface evolution or the rigid-lid:
36    C     |   U[n] = U* + dt x d/dx P
37    C     |   V[n] = V* + dt x d/dy P
38    C     |
39    C     | "Calculation of Gs"
40    C     | ===================
41    C     | This is where all the accelerations and tendencies (ie.
42    C     | physics, parameterizations etc...) are calculated
43    C     |   rho = rho ( theta[n], salt[n] )
44    C     |   b   = b(rho, theta)
45    C     |   K31 = K31 ( rho )
46    C     |   Gu[n] = Gu( u[n], v[n], wVel, b, ... )
47    C     |   Gv[n] = Gv( u[n], v[n], wVel, b, ... )
48    C     |   Gt[n] = Gt( theta[n], u[n], v[n], wVel, K31, ... )
49    C     |   Gs[n] = Gs( salt[n], u[n], v[n], wVel, K31, ... )
50    C     |
51    C     | "Time-stepping" or "Prediction"
52    C     | ================================
53    C     | The models variables are stepped forward with the appropriate
54    C     | time-stepping scheme (currently we use Adams-Bashforth II)
55    C     | - For momentum, the result is always *only* a "prediction"
56    C     | in that the flow may be divergent and will be "corrected"
57    C     | later with a surface pressure gradient.
58    C     | - Normally for tracers the result is the new field at time
59    C     | level [n+1} *BUT* in the case of implicit diffusion the result
60    C     | is also *only* a prediction.
61    C     | - We denote "predictors" with an asterisk (*).
62    C     |   U* = U[n] + dt x ( 3/2 Gu[n] - 1/2 Gu[n-1] )
63    C     |   V* = V[n] + dt x ( 3/2 Gv[n] - 1/2 Gv[n-1] )
64    C     |   theta[n+1] = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
65    C     |   salt[n+1] = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
66    C     | With implicit diffusion:
67    C     |   theta* = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
68    C     |   salt* = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )
69    C     |   (1 + dt * K * d_zz) theta[n] = theta*
70    C     |   (1 + dt * K * d_zz) salt[n] = salt*
71    C     |
72    C     *==========================================================*
73    C     \ev
74    
75    C     !USES:
76          IMPLICIT NONE
77  C     == Global variables ===  C     == Global variables ===
78  #include "SIZE.h"  #include "SIZE.h"
79  #include "EEPARAMS.h"  #include "EEPARAMS.h"
80  #include "PARAMS.h"  #include "PARAMS.h"
81  #include "DYNVARS.h"  #include "DYNVARS.h"
82  #include "GRID.h"  #include "GRID.h"
83    #include "GAD.h"
84  #ifdef ALLOW_PASSIVE_TRACER  #ifdef ALLOW_PASSIVE_TRACER
85  #include "TR1.h"  #include "TR1.h"
86  #endif  #endif
   
87  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
88  # include "tamc.h"  # include "tamc.h"
89  # include "tamc_keys.h"  # include "tamc_keys.h"
90  # include "FFIELDS.h"  # include "FFIELDS.h"
91    # include "EOS.h"
92  # ifdef ALLOW_KPP  # ifdef ALLOW_KPP
93  #  include "KPP.h"  #  include "KPP.h"
94  # endif  # endif
# Line 33  C     == Global variables === Line 96  C     == Global variables ===
96  #  include "GMREDI.h"  #  include "GMREDI.h"
97  # endif  # endif
98  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
   
99  #ifdef ALLOW_TIMEAVE  #ifdef ALLOW_TIMEAVE
100  #include "TIMEAVE_STATV.h"  #include "TIMEAVE_STATV.h"
101  #endif  #endif
102    
103    C     !INPUT/OUTPUT PARAMETERS:
104  C     == Routine arguments ==  C     == Routine arguments ==
105  C     myTime - Current time in simulation  C     myTime - Current time in simulation
106  C     myIter - Current iteration number in simulation  C     myIter - Current iteration number in simulation
# Line 46  C     myThid - Thread number for this in Line 109  C     myThid - Thread number for this in
109        INTEGER myIter        INTEGER myIter
110        INTEGER myThid        INTEGER myThid
111    
112    C     !LOCAL VARIABLES:
113  C     == Local variables  C     == Local variables
114  C     xA, yA                 - Per block temporaries holding face areas  C     xA, yA                 - Per block temporaries holding face areas
115  C     uTrans, vTrans, rTrans - Per block temporaries holding flow  C     uTrans, vTrans, rTrans - Per block temporaries holding flow
# Line 74  C     bi, bj Line 138  C     bi, bj
138  C     k, kup,        - Index for layer above and below. kup and kDown  C     k, kup,        - Index for layer above and below. kup and kDown
139  C     kDown, km1       are switched with layer to be the appropriate  C     kDown, km1       are switched with layer to be the appropriate
140  C                      index into fVerTerm.  C                      index into fVerTerm.
 C     tauAB - Adams-Bashforth timestepping weight: 0=forward ; 1/2=Adams-Bashf.  
141        _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
142        _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
143        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 84  C     tauAB - Adams-Bashforth timesteppi Line 147  C     tauAB - Adams-Bashforth timesteppi
147        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
148        _RL fVerS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
149        _RL fVerTr1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerTr1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
       _RL fVerU   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)  
       _RL fVerV   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)  
150        _RL phiHyd  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL phiHyd  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
151        _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
152        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 93  C     tauAB - Adams-Bashforth timesteppi Line 154  C     tauAB - Adams-Bashforth timesteppi
154        _RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiSurfY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
155        _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRT (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
156        _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL KappaRS (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
       _RL KappaRU (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)  
       _RL KappaRV (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)  
157        _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
158        _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
159        _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
160        _RL tauAB  C     This is currently used by IVDC and Diagnostics
   
 C This is currently used by IVDC and Diagnostics  
161        _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
   
162        INTEGER iMin, iMax        INTEGER iMin, iMax
163        INTEGER jMin, jMax        INTEGER jMin, jMax
164        INTEGER bi, bj        INTEGER bi, bj
165        INTEGER i, j        INTEGER i, j
166        INTEGER k, km1, kup, kDown        INTEGER k, km1, kup, kDown
167    
168  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)  
169    
 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---  
   
170  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
171  C--   dummy statement to end declaration part  C--   dummy statement to end declaration part
172        ikey = 1        ikey = 1
173          itdkey = 1
174  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
175    
176  C--   Set up work arrays with valid (i.e. not NaN) values  C--   Set up work arrays with valid (i.e. not NaN) values
# Line 175  C     uninitialised but inert locations. Line 184  C     uninitialised but inert locations.
184          yA(i,j)      = 0. _d 0          yA(i,j)      = 0. _d 0
185          uTrans(i,j)  = 0. _d 0          uTrans(i,j)  = 0. _d 0
186          vTrans(i,j)  = 0. _d 0          vTrans(i,j)  = 0. _d 0
         DO k=1,Nr  
          phiHyd(i,j,k)  = 0. _d 0  
          KappaRU(i,j,k) = 0. _d 0  
          KappaRV(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  
187          rhok   (i,j) = 0. _d 0          rhok   (i,j) = 0. _d 0
188          phiSurfX(i,j) = 0. _d 0          phiSurfX(i,j) = 0. _d 0
189          phiSurfY(i,j) = 0. _d 0          phiSurfY(i,j) = 0. _d 0
# Line 200  CHPF$ INDEPENDENT Line 200  CHPF$ INDEPENDENT
200    
201  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
202  C--    HPF directive to help TAMC  C--    HPF directive to help TAMC
203  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS,fVerU,fVerV  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS
204  CHPF$&                  ,phiHyd,utrans,vtrans,xA,yA  CHPF$&                  ,phiHyd,utrans,vtrans,xA,yA
205  CHPF$&                  ,KappaRT,KappaRS,KappaRU,KappaRV  CHPF$&                  ,KappaRT,KappaRS
206  CHPF$&                  )  CHPF$&                  )
207  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
208    
# Line 211  CHPF$&                  ) Line 211  CHPF$&                  )
211  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
212            act1 = bi - myBxLo(myThid)            act1 = bi - myBxLo(myThid)
213            max1 = myBxHi(myThid) - myBxLo(myThid) + 1            max1 = myBxHi(myThid) - myBxLo(myThid) + 1
   
214            act2 = bj - myByLo(myThid)            act2 = bj - myByLo(myThid)
215            max2 = myByHi(myThid) - myByLo(myThid) + 1            max2 = myByHi(myThid) - myByLo(myThid) + 1
   
216            act3 = myThid - 1            act3 = myThid - 1
217            max3 = nTx*nTy            max3 = nTx*nTy
   
218            act4 = ikey_dynamics - 1            act4 = ikey_dynamics - 1
219              itdkey = (act1 + 1) + act2*max1
           ikey = (act1 + 1) + act2*max1  
220       &                      + act3*max1*max2       &                      + act3*max1*max2
221       &                      + act4*max1*max2*max3       &                      + act4*max1*max2*max3
222  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
# Line 235  C--     Set up work arrays that need val Line 231  C--     Set up work arrays that need val
231            fVerS  (i,j,2) = 0. _d 0            fVerS  (i,j,2) = 0. _d 0
232            fVerTr1(i,j,1) = 0. _d 0            fVerTr1(i,j,1) = 0. _d 0
233            fVerTr1(i,j,2) = 0. _d 0            fVerTr1(i,j,2) = 0. _d 0
234            fVerU  (i,j,1) = 0. _d 0            rhoKM1 (i,j)   = 0. _d 0
           fVerU  (i,j,2) = 0. _d 0  
           fVerV  (i,j,1) = 0. _d 0  
           fVerV  (i,j,2) = 0. _d 0  
235           ENDDO           ENDDO
236          ENDDO          ENDDO
237    
# Line 246  C--     Set up work arrays that need val Line 239  C--     Set up work arrays that need val
239           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
240            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
241  C This is currently also used by IVDC and Diagnostics  C This is currently also used by IVDC and Diagnostics
242               phiHyd(i,j,k) = 0. _d 0
243               sigmaX(i,j,k) = 0. _d 0
244               sigmaY(i,j,k) = 0. _d 0
245               sigmaR(i,j,k) = 0. _d 0
246             ConvectCount(i,j,k) = 0.             ConvectCount(i,j,k) = 0.
247             KappaRT(i,j,k) = 0. _d 0             KappaRT(i,j,k)    = 0. _d 0
248             KappaRS(i,j,k) = 0. _d 0             KappaRS(i,j,k)    = 0. _d 0
249    #ifdef ALLOW_AUTODIFF_TAMC
250    cph all the following init. are necessary for TAF
251    cph although some of these are re-initialised later.
252               gT(i,j,k,bi,bj)   = 0. _d 0
253               gS(i,j,k,bi,bj)   = 0. _d 0
254    # ifdef ALLOW_PASSIVE_TRACER
255               gTr1(i,j,k,bi,bj) = 0. _d 0
256    # endif
257    # ifdef ALLOW_GMREDI
258               Kwx(i,j,k,bi,bj)  = 0. _d 0
259               Kwy(i,j,k,bi,bj)  = 0. _d 0
260               Kwz(i,j,k,bi,bj)  = 0. _d 0
261    #  ifdef GM_NON_UNITY_DIAGONAL
262               Kux(i,j,k,bi,bj)  = 0. _d 0
263               Kvy(i,j,k,bi,bj)  = 0. _d 0
264    #  endif
265    #  ifdef GM_EXTRA_DIAGONAL
266               Kuz(i,j,k,bi,bj)  = 0. _d 0
267               Kvz(i,j,k,bi,bj)  = 0. _d 0
268    #  endif
269    #  ifdef GM_BOLUS_ADVEC
270               GM_PsiX(i,j,k,bi,bj)  = 0. _d 0
271               GM_PsiY(i,j,k,bi,bj)  = 0. _d 0
272    #  endif
273    # endif /* ALLOW_GMREDI */
274    #endif /* ALLOW_AUTODIFF_TAMC */
275            ENDDO            ENDDO
276           ENDDO           ENDDO
277          ENDDO          ENDDO
278    
279          iMin = 1-OLx+1          iMin = 1-OLx
280          iMax = sNx+OLx          iMax = sNx+OLx
281          jMin = 1-OLy+1          jMin = 1-OLy
282          jMax = sNy+OLy          jMax = sNy+OLy
283    
284    
285  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
286  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
287  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
288  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  #ifdef ALLOW_KPP
289  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
290  #ifdef ALLOW_PASSIVE_TRACER  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
 CADJ STORE tr1  (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  
291  #endif  #endif
292  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
293    
# Line 277  C? Patrick, is this formula correct now Line 299  C? Patrick, is this formula correct now
299  C? Do we still need this?  C? Do we still need this?
300  cph kkey formula corrected.  cph kkey formula corrected.
301  cph Needed for rhok, rhokm1, in the case useGMREDI.  cph Needed for rhok, rhokm1, in the case useGMREDI.
302           kkey = (ikey-1)*Nr + k           kkey = (itdkey-1)*Nr + k
 CADJ STORE rhokm1(:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte  
 CADJ STORE rhok  (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte  
303  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
304    
305  C--       Integrate continuity vertically for vertical velocity  C--       Integrate continuity vertically for vertical velocity
306            CALL INTEGRATE_FOR_W(  c         CALL INTEGRATE_FOR_W(
307       I                         bi, bj, k, uVel, vVel,  c    I                         bi, bj, k, uVel, vVel,
308       O                         wVel,  c    O                         wVel,
309       I                         myThid )  c    I                         myThid )
310    
311  #ifdef    ALLOW_OBCS  #ifdef    ALLOW_OBCS
312  #ifdef    ALLOW_NONHYDROSTATIC  #ifdef    ALLOW_NONHYDROSTATIC
313  C--       Apply OBC to W if in N-H mode  C--       Apply OBC to W if in N-H mode
314            IF (useOBCS.AND.nonHydrostatic) THEN  c         IF (useOBCS.AND.nonHydrostatic) THEN
315              CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )  c           CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )
316            ENDIF  c         ENDIF
317  #endif    /* ALLOW_NONHYDROSTATIC */  #endif    /* ALLOW_NONHYDROSTATIC */
318  #endif    /* ALLOW_OBCS */  #endif    /* ALLOW_OBCS */
319    
320    C--       Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h
321    C--       MOST of THERMODYNAMICS will be disabled
322    #ifndef SINGLE_LAYER_MODE
323    
324  C--       Calculate gradients of potential density for isoneutral  C--       Calculate gradients of potential density for isoneutral
325  C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)  C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)
326  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 304  c         IF ( k.GT.1 .AND. (useGMRedi.O Line 328  c         IF ( k.GT.1 .AND. (useGMRedi.O
328  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
329  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
330  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
331    CADJ STORE pressure(:,:,k,bi,bj) =
332    CADJ &     comlev1_bibj_k, key=kkey, byte=isbyte
333  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
334              CALL FIND_RHO(              CALL FIND_RHO(
335       I        bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,       I        bi, bj, iMin, iMax, jMin, jMax, k, k,
336       I        theta, salt,       I        theta, salt,
337       O        rhoK,       O        rhoK,
338       I        myThid )       I        myThid )
339    
340              IF (k.GT.1) THEN              IF (k.GT.1) THEN
341  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
342  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
343  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
344    CADJ STORE pressure(:,:,k-1,bi,bj) =
345    CADJ &     comlev1_bibj_k, key=kkey, byte=isbyte
346  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
347               CALL FIND_RHO(               CALL FIND_RHO(
348       I        bi, bj, iMin, iMax, jMin, jMax, k-1, k, eosType,       I        bi, bj, iMin, iMax, jMin, jMax, k-1, k,
349       I        theta, salt,       I        theta, salt,
350       O        rhoKm1,       O        rhoKm1,
351       I        myThid )       I        myThid )
# Line 328  CADJ STORE salt (:,:,k-1,bi,bj) = comlev Line 357  CADJ STORE salt (:,:,k-1,bi,bj) = comlev
357       I             myThid )       I             myThid )
358            ENDIF            ENDIF
359    
360    #ifdef ALLOW_AUTODIFF_TAMC
361    CADJ STORE rhok   (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte
362    CADJ STORE rhokm1 (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte
363    #endif /* ALLOW_AUTODIFF_TAMC */
364  C--       Implicit Vertical Diffusion for Convection  C--       Implicit Vertical Diffusion for Convection
365  c ==> should use sigmaR !!!  c ==> should use sigmaR !!!
366            IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN            IF (k.GT.1 .AND. ivdc_kappa.NE.0.) THEN
# Line 338  c ==> should use sigmaR !!! Line 371  c ==> should use sigmaR !!!
371       I        myTime, myIter, myThid)       I        myTime, myIter, myThid)
372            ENDIF            ENDIF
373    
374    #endif /* SINGLE_LAYER_MODE */
375    
376  C--     end of diagnostic k loop (Nr:1)  C--     end of diagnostic k loop (Nr:1)
377          ENDDO          ENDDO
378    
379  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
380  cph avoids recomputation of integrate_for_w  cph avoids recomputation of integrate_for_w
381  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
382    CADJ STORE pressure (:,:,:,bi,bj) =
383    CADJ &     comlev1_bibj, key=itdkey, byte=isbyte
384  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
385    
386  #ifdef  ALLOW_OBCS  #ifdef  ALLOW_OBCS
387  C--     Calculate future values on open boundaries  C--     Calculate future values on open boundaries
388          IF (useOBCS) THEN          IF (useOBCS) THEN
389            CALL OBCS_CALC( bi, bj, myTime+deltaT,            CALL OBCS_CALC( bi, bj, myTime+deltaT, myIter+1,
390       I            uVel, vVel, wVel, theta, salt,       I            uVel, vVel, wVel, theta, salt,
391       I            myThid )       I            myThid )
392          ENDIF          ENDIF
393  #endif  /* ALLOW_OBCS */  #endif  /* ALLOW_OBCS */
394    
395    
396    c********************************************
397    cswdice --- add ---
398    #ifdef ALLOW_THERM_SEAICE
399    C--     Determines forcing terms based on external fields
400    c--     including effects from ice
401            CALL ICE_FORCING(
402         I             bi, bj, iMin, iMax, jMin, jMax,
403         I             myThid )
404    #else
405    
406    cswdice --- end add ---
407    
408  C--     Determines forcing terms based on external fields  C--     Determines forcing terms based on external fields
409  C       relaxation terms, etc.  C       relaxation terms, etc.
410          CALL EXTERNAL_FORCING_SURF(          CALL EXTERNAL_FORCING_SURF(
411       I             bi, bj, iMin, iMax, jMin, jMax,       I             bi, bj, iMin, iMax, jMin, jMax,
412       I             myThid )       I             myThid )
413    cswdice --- add ----
414    #endif
415    cswdice --- end add ---
416    c******************************************
417    
418    
419    
420    
421  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
422  cph needed for KPP  cph needed for KPP
423  CADJ STORE surfacetendencyU(:,:,bi,bj)  CADJ STORE surfacetendencyU(:,:,bi,bj)
424  CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
425  CADJ STORE surfacetendencyV(:,:,bi,bj)  CADJ STORE surfacetendencyV(:,:,bi,bj)
426  CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
427  CADJ STORE surfacetendencyS(:,:,bi,bj)  CADJ STORE surfacetendencyS(:,:,bi,bj)
428  CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
429  CADJ STORE surfacetendencyT(:,:,bi,bj)  CADJ STORE surfacetendencyT(:,:,bi,bj)
430  CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
431  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
432    
433    C--     Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h
434    C--     MOST of THERMODYNAMICS will be disabled
435    #ifndef SINGLE_LAYER_MODE
436    
437  #ifdef  ALLOW_GMREDI  #ifdef  ALLOW_GMREDI
438    
439  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
440  CADJ STORE sigmaX(:,:,:) = comlev1, key=ikey, byte=isbyte  cph storing here is needed only for one GMREDI_OPTIONS:
441  CADJ STORE sigmaY(:,:,:) = comlev1, key=ikey, byte=isbyte  cph define GM_BOLUS_ADVEC
442  CADJ STORE sigmaR(:,:,:) = comlev1, key=ikey, byte=isbyte  cph but I've avoided the #ifdef for now, in case more things change
443    CADJ STORE sigmaX(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte
444    CADJ STORE sigmaY(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte
445    CADJ STORE sigmaR(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte
446  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
447  C--     Calculate iso-neutral slopes for the GM/Redi parameterisation  C--     Calculate iso-neutral slopes for the GM/Redi parameterisation
448          IF (useGMRedi) THEN          IF (useGMRedi) THEN
449            DO k=1,Nr            CALL GMREDI_CALC_TENSOR(
450              CALL GMREDI_CALC_TENSOR(       I             bi, bj, iMin, iMax, jMin, jMax,
      I             bi, bj, iMin, iMax, jMin, jMax, k,  
451       I             sigmaX, sigmaY, sigmaR,       I             sigmaX, sigmaY, sigmaR,
452       I             myThid )       I             myThid )
           ENDDO  
453  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
454          ELSE          ELSE
455            DO k=1, Nr            CALL GMREDI_CALC_TENSOR_DUMMY(
456              CALL GMREDI_CALC_TENSOR_DUMMY(       I             bi, bj, iMin, iMax, jMin, jMax,
      I             bi, bj, iMin, iMax, jMin, jMax, k,  
457       I             sigmaX, sigmaY, sigmaR,       I             sigmaX, sigmaY, sigmaR,
458       I             myThid )       I             myThid )
           ENDDO  
459  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
460          ENDIF          ENDIF
461    
462  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
463  CADJ STORE Kwx(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE Kwx(:,:,:,bi,bj)   = comlev1_bibj, key=itdkey, byte=isbyte
464  CADJ STORE Kwy(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE Kwy(:,:,:,bi,bj)   = comlev1_bibj, key=itdkey, byte=isbyte
465  CADJ STORE Kwz(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE Kwz(:,:,:,bi,bj)   = comlev1_bibj, key=itdkey, byte=isbyte
466  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
467    
468  #endif  /* ALLOW_GMREDI */  #endif  /* ALLOW_GMREDI */
# Line 420  C--     Compute KPP mixing coefficients Line 481  C--     Compute KPP mixing coefficients
481    
482  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
483  CADJ STORE KPPghat   (:,:,:,bi,bj)  CADJ STORE KPPghat   (:,:,:,bi,bj)
 CADJ &   , KPPviscAz (:,:,:,bi,bj)  
484  CADJ &   , KPPdiffKzT(:,:,:,bi,bj)  CADJ &   , KPPdiffKzT(:,:,:,bi,bj)
485  CADJ &   , KPPdiffKzS(:,:,:,bi,bj)  CADJ &   , KPPdiffKzS(:,:,:,bi,bj)
486  CADJ &   , KPPfrac   (:,:  ,bi,bj)  CADJ &   , KPPfrac   (:,:  ,bi,bj)
487  CADJ &                 = comlev1_bibj, key=ikey, byte=isbyte  CADJ &                 = comlev1_bibj, key=itdkey, byte=isbyte
488  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
489    
490  #endif  /* ALLOW_KPP */  #endif  /* ALLOW_KPP */
491    
492  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
493  CADJ STORE KappaRT(:,:,:)     = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE KappaRT(:,:,:)     = comlev1_bibj, key=itdkey, byte=isbyte
494  CADJ STORE KappaRS(:,:,:)     = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE KappaRS(:,:,:)     = comlev1_bibj, key=itdkey, byte=isbyte
495  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
496  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
497  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
498  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
499  #ifdef ALLOW_PASSIVE_TRACER  #ifdef ALLOW_PASSIVE_TRACER
500  CADJ STORE tr1  (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE tr1  (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
501  #endif  #endif
502  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
503    
504  #ifdef ALLOW_AIM  #ifdef ALLOW_AIM
505  C       AIM - atmospheric intermediate model, physics package code.  C       AIM - atmospheric intermediate model, physics package code.
 C note(jmc) : phiHyd=0 at this point but is not really used in Molteni Physics  
506          IF ( useAIM ) THEN          IF ( useAIM ) THEN
507           CALL TIMER_START('AIM_DO_ATMOS_PHYS      [DYNAMICS]', myThid)           CALL TIMER_START('AIM_DO_PHYSICS   [THERMODYNAMICS]', myThid)
508           CALL AIM_DO_ATMOS_PHYSICS( phiHyd, bi, bj, myTime, myThid )           CALL AIM_DO_PHYSICS( bi, bj, myTime, myIter, myThid )
509           CALL TIMER_STOP ('AIM_DO_ATMOS_PHYS      [DYNAMICS]', myThid)           CALL TIMER_STOP( 'AIM_DO_PHYSICS   [THERMODYNAMICS]', myThid)
510          ENDIF          ENDIF
511  #endif /* ALLOW_AIM */  #endif /* ALLOW_AIM */
512    
513    #ifdef ALLOW_TIMEAVE
514            IF (taveFreq.GT.0. .AND. ivdc_kappa.NE.0.) THEN
515              CALL TIMEAVE_CUMULATE(ConvectCountTave, ConvectCount, Nr,
516         I                               deltaTclock, bi, bj, myThid)
517            ENDIF
518    #endif /* ALLOW_TIMEAVE */
519    
520    #ifndef DISABLE_MULTIDIM_ADVECTION
521    C--     Some advection schemes are better calculated using a multi-dimensional
522    C       method in the absence of any other terms and, if used, is done here.
523    C
524    C The CPP flag DISABLE_MULTIDIM_ADVECTION is currently unset in GAD_OPTIONS.h
525    C The default is to use multi-dimensinal advection for non-linear advection
526    C schemes. However, for the sake of efficiency of the adjoint it is necessary
527    C to be able to exclude this scheme to avoid excessive storage and
528    C recomputation. It *is* differentiable, if you need it.
529    C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to
530    C disable this section of code.
531            IF (tempMultiDimAdvec) THEN
532              CALL GAD_ADVECTION(bi,bj,tempAdvScheme,GAD_TEMPERATURE,
533         U                      theta,gT,
534         I                      myTime,myIter,myThid)
535            ENDIF
536            IF (saltMultiDimAdvec) THEN
537              CALL GAD_ADVECTION(bi,bj,saltAdvScheme,GAD_SALINITY,
538         U                      salt,gS,
539         I                      myTime,myIter,myThid)
540            ENDIF
541    C Since passive tracers are configurable separately from T,S we
542    C call the multi-dimensional method for PTRACERS regardless
543    C of whether multiDimAdvection is set or not.
544    #ifdef ALLOW_PTRACERS
545            IF ( usePTRACERS ) THEN
546             CALL PTRACERS_ADVECTION( bi,bj,myIter,myTime,myThid )
547            ENDIF
548    #endif /* ALLOW_PTRACERS */
549    #endif /* DISABLE_MULTIDIM_ADVECTION */
550    
551  C--     Start of thermodynamics loop  C--     Start of thermodynamics loop
552          DO k=Nr,1,-1          DO k=Nr,1,-1
# Line 458  C--     Start of thermodynamics loop Line 554  C--     Start of thermodynamics loop
554  C? Patrick Is this formula correct?  C? Patrick Is this formula correct?
555  cph Yes, but I rewrote it.  cph Yes, but I rewrote it.
556  cph Also, the KappaR? need the index and subscript k!  cph Also, the KappaR? need the index and subscript k!
557           kkey = (ikey-1)*Nr + k           kkey = (itdkey-1)*Nr + k
558  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
559    
560  C--       km1    Points to level above k (=k-1)  C--       km1    Points to level above k (=k-1)
# Line 480  C--       Get temporary terms used by te Line 576  C--       Get temporary terms used by te
576       O         xA,yA,uTrans,vTrans,rTrans,maskUp,       O         xA,yA,uTrans,vTrans,rTrans,maskUp,
577       I         myThid)       I         myThid)
578    
579    #ifdef ALLOW_GMREDI
580    C--   Residual transp = Bolus transp + Eulerian transp
581              IF (useGMRedi) THEN
582                CALL GMREDI_CALC_UVFLOW(
583         &            uTrans, vTrans, bi, bj, k, myThid)
584                IF (K.GE.2) CALL GMREDI_CALC_WFLOW(
585         &                    rTrans, bi, bj, k, myThid)
586              ENDIF
587    #endif /* ALLOW_GMREDI */
588    
589  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
590  CADJ STORE KappaRT(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE KappaRT(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte
591  CADJ STORE KappaRS(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE KappaRS(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte
# Line 490  C--      Calculate the total vertical di Line 596  C--      Calculate the total vertical di
596           CALL CALC_DIFFUSIVITY(           CALL CALC_DIFFUSIVITY(
597       I        bi,bj,iMin,iMax,jMin,jMax,k,       I        bi,bj,iMin,iMax,jMin,jMax,k,
598       I        maskUp,       I        maskUp,
599       O        KappaRT,KappaRS,KappaRU,KappaRV,       O        KappaRT,KappaRS,
600       I        myThid)       I        myThid)
601  #endif  #endif
602    
# Line 507  C        and step forward storing result Line 613  C        and step forward storing result
613       I         xA,yA,uTrans,vTrans,rTrans,maskUp,       I         xA,yA,uTrans,vTrans,rTrans,maskUp,
614       I         KappaRT,       I         KappaRT,
615       U         fVerT,       U         fVerT,
616       I         myTime, myThid)       I         myTime,myIter,myThid)
            tauAB = 0.5d0 + abEps  
617             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
618       I         bi,bj,iMin,iMax,jMin,jMax,k,tauAB,       I         bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,
619       I         theta, gT,       I         theta, gT,
      U         gTnm1,  
620       I         myIter, myThid)       I         myIter, myThid)
621           ENDIF           ENDIF
622    cswdice ---- add ---
623    #ifdef ALLOW_THERM_SEAICE
624           if (k.eq.1) then
625            call ICE_FREEZE(bi, bj, iMin, iMax, jMin, jMax, myThid )
626           endif
627    #endif
628    cswdice -- end add ---
629           IF ( saltStepping ) THEN           IF ( saltStepping ) THEN
630             CALL CALC_GS(             CALL CALC_GS(
631       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
632       I         xA,yA,uTrans,vTrans,rTrans,maskUp,       I         xA,yA,uTrans,vTrans,rTrans,maskUp,
633       I         KappaRS,       I         KappaRS,
634       U         fVerS,       U         fVerS,
635       I         myTime, myThid)       I         myTime,myIter,myThid)
            tauAB = 0.5d0 + abEps  
636             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
637       I         bi,bj,iMin,iMax,jMin,jMax,k,tauAB,       I         bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,
638       I         salt, gS,       I         salt, gS,
      U         gSnm1,  
639       I         myIter, myThid)       I         myIter, myThid)
640           ENDIF           ENDIF
641  #ifdef ALLOW_PASSIVE_TRACER  #ifdef ALLOW_PASSIVE_TRACER
# Line 536  C        and step forward storing result Line 645  C        and step forward storing result
645       I         xA,yA,uTrans,vTrans,rTrans,maskUp,       I         xA,yA,uTrans,vTrans,rTrans,maskUp,
646       I         KappaRT,       I         KappaRT,
647       U         fVerTr1,       U         fVerTr1,
648       I         myTime, myThid)       I         myTime,myIter,myThid)
            tauAB = 0.5d0 + abEps  
649             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
650       I         bi,bj,iMin,iMax,jMin,jMax,k,tauAB,       I         bi,bj,iMin,iMax,jMin,jMax,k,tracerAdvScheme,
651       I         Tr1, gTr1,       I         Tr1, gTr1,
652       U         gTr1NM1,       I         myIter,myThid)
      I         myIter, myThid)  
653           ENDIF           ENDIF
654  #endif  #endif
655    #ifdef ALLOW_PTRACERS
656             IF ( usePTRACERS ) THEN
657               CALL PTRACERS_INTEGERATE(
658         I         bi,bj,k,
659         I         xA,yA,uTrans,vTrans,rTrans,maskUp,
660         X         KappaRS,
661         I         myIter,myTime,myThid)
662             ENDIF
663    #endif /* ALLOW_PTRACERS */
664    
665  #ifdef   ALLOW_OBCS  #ifdef   ALLOW_OBCS
666  C--      Apply open boundary conditions  C--      Apply open boundary conditions
667           IF (useOBCS) THEN           IF (useOBCS) THEN
668             CALL OBCS_APPLY_TS( bi, bj, k, gTnm1, gSnm1, myThid )             CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )
669           END IF           END IF
670  #endif   /* ALLOW_OBCS */  #endif   /* ALLOW_OBCS */
671    
672  C--      Freeze water  C--      Freeze water
673           IF (allowFreezing) THEN           IF ( allowFreezing .AND. .NOT. useSEAICE ) THEN
674  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
675  CADJ STORE gTNm1(:,:,k,bi,bj) = comlev1_bibj_k  CADJ STORE gT(:,:,k,bi,bj) = comlev1_bibj_k
676  CADJ &   , key = kkey, byte = isbyte  CADJ &   , key = kkey, byte = isbyte
677  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
678              CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )              CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )
# Line 565  CADJ &   , key = kkey, byte = isbyte Line 681  CADJ &   , key = kkey, byte = isbyte
681  C--     end of thermodynamic k loop (Nr:1)  C--     end of thermodynamic k loop (Nr:1)
682          ENDDO          ENDDO
683    
684    cswdice -- add ---
685    #ifdef ALLOW_THERM_SEAICE
686    c timeaveraging for ice model values
687               CALL ICE_AVE(bi,bj,iMin,iMax,jMin,jMax,myThid )
688    #endif
689    cswdice --- end add ---
690    
691    
692    
 #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 */  
693    
694  C--     Implicit diffusion  C--     Implicit diffusion
695          IF (implicitDiffusion) THEN          IF (implicitDiffusion) THEN
696    
697           IF (tempStepping) THEN           IF (tempStepping) THEN
698  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
699              idkey = iikey + 1  CADJ STORE gT(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
 CADJ STORE gTNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  
700  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
701              CALL IMPLDIFF(              CALL IMPLDIFF(
702       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
703       I         deltaTtracer, KappaRT, recip_HFacC,       I         deltaTtracer, KappaRT, recip_HFacC,
704       U         gTNm1,       U         gT,
705       I         myThid )       I         myThid )
706           ENDIF           ENDIF
707    
708           IF (saltStepping) THEN           IF (saltStepping) THEN
709  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
710           idkey = iikey + 2  CADJ STORE gS(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
 CADJ STORE gSNm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  
711  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
712              CALL IMPLDIFF(              CALL IMPLDIFF(
713       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
714       I         deltaTtracer, KappaRS, recip_HFacC,       I         deltaTtracer, KappaRS, recip_HFacC,
715       U         gSNm1,       U         gS,
716       I         myThid )       I         myThid )
717           ENDIF           ENDIF
718    
719  #ifdef ALLOW_PASSIVE_TRACER  #ifdef ALLOW_PASSIVE_TRACER
720           IF (tr1Stepping) THEN           IF (tr1Stepping) THEN
721  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
722  CADJ STORE gTr1Nm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE gTr1(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
723  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
724            CALL IMPLDIFF(            CALL IMPLDIFF(
725       I      bi, bj, iMin, iMax, jMin, jMax,       I      bi, bj, iMin, iMax, jMin, jMax,
726       I      deltaTtracer, KappaRT, recip_HFacC,       I      deltaTtracer, KappaRT, recip_HFacC,
727       U      gTr1Nm1,       U      gTr1,
728       I      myThid )       I      myThid )
729           ENDIF           ENDIF
730  #endif  #endif
731    
732    #ifdef ALLOW_PTRACERS
733    C Vertical diffusion (implicit) for passive tracers
734             IF ( usePTRACERS ) THEN
735               CALL PTRACERS_IMPLDIFF( bi,bj,KappaRS,myThid )
736             ENDIF
737    #endif /* ALLOW_PTRACERS */
738    
739  #ifdef   ALLOW_OBCS  #ifdef   ALLOW_OBCS
740  C--      Apply open boundary conditions  C--      Apply open boundary conditions
741           IF (useOBCS) THEN           IF (useOBCS) THEN
742             DO K=1,Nr             DO K=1,Nr
743               CALL OBCS_APPLY_TS( bi, bj, k, gTnm1, gSnm1, myThid )               CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )
744             ENDDO             ENDDO
745           END IF           END IF
746  #endif   /* ALLOW_OBCS */  #endif   /* ALLOW_OBCS */
# Line 628  C--      Apply open boundary conditions Line 748  C--      Apply open boundary conditions
748  C--     End If implicitDiffusion  C--     End If implicitDiffusion
749          ENDIF          ENDIF
750    
751    #endif /* SINGLE_LAYER_MODE */
752    
753  Ccs-  Ccs-
754         ENDDO         ENDDO
755        ENDDO        ENDDO
756    
757  #ifdef ALLOW_AIM  #ifdef ALLOW_AIM
758        IF ( useAIM ) THEN  c     IF ( useAIM ) THEN
759         CALL AIM_AIM2DYN_EXCHANGES( myTime, myThid )  c      CALL AIM_AIM2DYN_EXCHANGES( myTime, myThid )
760        ENDIF  c     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  
761  #endif /* ALLOW_AIM */  #endif /* ALLOW_AIM */
762    c     IF ( staggerTimeStep ) THEN
763    c      IF ( useAIM .OR. useCubedSphereExchange ) THEN
764    c        IF (tempStepping) _EXCH_XYZ_R8(gT,myThid)
765    c        IF (saltStepping) _EXCH_XYZ_R8(gS,myThid)
766    c      ELSEIF ( useGMRedi .AND. Oly.LT.4 ) THEN
767    cc       .AND. GM_AdvForm .AND. .NOT.GM_AdvSeparate ) THEN
768    c        IF (tempMultiDimAdvec) _EXCH_XYZ_R8(gT,myThid)
769    c        IF (saltMultiDimAdvec) _EXCH_XYZ_R8(gS,myThid)
770    c      ENDIF  
771    c     ENDIF  
772    
773    #ifndef DISABLE_DEBUGMODE
774          If (debugMode) THEN
775           CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (THERMODYNAMICS)',myThid)
776           CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (THERMODYNAMICS)',myThid)
777           CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (THERMODYNAMICS)',myThid)
778           CALL DEBUG_STATS_RL(Nr,theta,'Theta (THERMODYNAMICS)',myThid)
779           CALL DEBUG_STATS_RL(Nr,salt,'Salt (THERMODYNAMICS)',myThid)
780           CALL DEBUG_STATS_RL(Nr,Gt,'Gt (THERMODYNAMICS)',myThid)
781           CALL DEBUG_STATS_RL(Nr,Gs,'Gs (THERMODYNAMICS)',myThid)
782           CALL DEBUG_STATS_RL(Nr,GtNm1,'GtNm1 (THERMODYNAMICS)',myThid)
783           CALL DEBUG_STATS_RL(Nr,GsNm1,'GsNm1 (THERMODYNAMICS)',myThid)
784    #ifdef ALLOW_PTRACERS
785           IF ( usePTRACERS ) THEN
786             CALL PTRACERS_DEBUG(myThid)
787           ENDIF
788    #endif /* ALLOW_PTRACERS */
789          ENDIF
790    #endif
791    
792        RETURN        RETURN
793        END        END
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.34
  ViewVC Help
Powered by ViewVC 1.1.22