/[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.2 by heimbach, Mon Aug 13 18:05:26 2001 UTC revision 1.35 by heimbach, Fri Jan 10 23:41:15 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 94  C     tauAB - Adams-Bashforth timesteppi Line 157  C     tauAB - Adams-Bashforth timesteppi
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 171  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  
          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 205  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 229  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              rhoKM1 (i,j)   = 0. _d 0
235           ENDDO           ENDDO
236          ENDDO          ENDDO
237    
# Line 236  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    #ifdef ALLOW_KPP
289    CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
290    CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
291    #endif
292  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
293    
294  C--     Start of diagnostic loop  C--     Start of diagnostic loop
# Line 262  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 289  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 313  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 323  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    
448  C--     Calculate iso-neutral slopes for the GM/Redi parameterisation  C--     Calculate iso-neutral slopes for the GM/Redi parameterisation
449          IF (useGMRedi) THEN          IF (useGMRedi) THEN
450            DO k=1,Nr            CALL GMREDI_CALC_TENSOR(
451              CALL GMREDI_CALC_TENSOR(       I             bi, bj, iMin, iMax, jMin, jMax,
      I             bi, bj, iMin, iMax, jMin, jMax, k,  
452       I             sigmaX, sigmaY, sigmaR,       I             sigmaX, sigmaY, sigmaR,
453       I             myThid )       I             myThid )
           ENDDO  
454  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
455          ELSE          ELSE
456            DO k=1, Nr            CALL GMREDI_CALC_TENSOR_DUMMY(
457              CALL GMREDI_CALC_TENSOR_DUMMY(       I             bi, bj, iMin, iMax, jMin, jMax,
      I             bi, bj, iMin, iMax, jMin, jMax, k,  
458       I             sigmaX, sigmaY, sigmaR,       I             sigmaX, sigmaY, sigmaR,
459       I             myThid )       I             myThid )
           ENDDO  
460  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
461          ENDIF          ENDIF
462    
463  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
464  CADJ STORE Kwx(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE Kwx(:,:,:,bi,bj)   = comlev1_bibj, key=itdkey, byte=isbyte
465  CADJ STORE Kwy(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE Kwy(:,:,:,bi,bj)   = comlev1_bibj, key=itdkey, byte=isbyte
466  CADJ STORE Kwz(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE Kwz(:,:,:,bi,bj)   = comlev1_bibj, key=itdkey, byte=isbyte
467  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
468    
469  #endif  /* ALLOW_GMREDI */  #endif  /* ALLOW_GMREDI */
# Line 405  C--     Compute KPP mixing coefficients Line 482  C--     Compute KPP mixing coefficients
482    
483  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
484  CADJ STORE KPPghat   (:,:,:,bi,bj)  CADJ STORE KPPghat   (:,:,:,bi,bj)
 CADJ &   , KPPviscAz (:,:,:,bi,bj)  
485  CADJ &   , KPPdiffKzT(:,:,:,bi,bj)  CADJ &   , KPPdiffKzT(:,:,:,bi,bj)
486  CADJ &   , KPPdiffKzS(:,:,:,bi,bj)  CADJ &   , KPPdiffKzS(:,:,:,bi,bj)
487  CADJ &   , KPPfrac   (:,:  ,bi,bj)  CADJ &   , KPPfrac   (:,:  ,bi,bj)
488  CADJ &                 = comlev1_bibj, key=ikey, byte=isbyte  CADJ &                 = comlev1_bibj, key=itdkey, byte=isbyte
489  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
490    
491  #endif  /* ALLOW_KPP */  #endif  /* ALLOW_KPP */
492    
493  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
494  CADJ STORE KappaRT(:,:,:)     = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE KappaRT(:,:,:)     = comlev1_bibj, key=itdkey, byte=isbyte
495  CADJ STORE KappaRS(:,:,:)     = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE KappaRS(:,:,:)     = comlev1_bibj, key=itdkey, byte=isbyte
496  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
497  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
498  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
499  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
500  #ifdef ALLOW_PASSIVE_TRACER  #ifdef ALLOW_PASSIVE_TRACER
501  CADJ STORE tr1  (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE tr1  (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
502  #endif  #endif
503  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
504    
505  #ifdef ALLOW_AIM  #ifdef ALLOW_AIM
506  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  
507          IF ( useAIM ) THEN          IF ( useAIM ) THEN
508           CALL TIMER_START('AIM_DO_ATMOS_PHYS      [DYNAMICS]', myThid)           CALL TIMER_START('AIM_DO_PHYSICS   [THERMODYNAMICS]', myThid)
509           CALL AIM_DO_ATMOS_PHYSICS( phiHyd, bi, bj, myTime, myThid )           CALL AIM_DO_PHYSICS( bi, bj, myTime, myIter, myThid )
510           CALL TIMER_STOP ('AIM_DO_ATMOS_PHYS      [DYNAMICS]', myThid)           CALL TIMER_STOP( 'AIM_DO_PHYSICS   [THERMODYNAMICS]', myThid)
511          ENDIF          ENDIF
512  #endif /* ALLOW_AIM */  #endif /* ALLOW_AIM */
513    
514    #ifdef ALLOW_TIMEAVE
515            IF (taveFreq.GT.0. .AND. ivdc_kappa.NE.0.) THEN
516              CALL TIMEAVE_CUMULATE(ConvectCountTave, ConvectCount, Nr,
517         I                               deltaTclock, bi, bj, myThid)
518            ENDIF
519    #endif /* ALLOW_TIMEAVE */
520    
521    #ifndef DISABLE_MULTIDIM_ADVECTION
522    C--     Some advection schemes are better calculated using a multi-dimensional
523    C       method in the absence of any other terms and, if used, is done here.
524    C
525    C The CPP flag DISABLE_MULTIDIM_ADVECTION is currently unset in GAD_OPTIONS.h
526    C The default is to use multi-dimensinal advection for non-linear advection
527    C schemes. However, for the sake of efficiency of the adjoint it is necessary
528    C to be able to exclude this scheme to avoid excessive storage and
529    C recomputation. It *is* differentiable, if you need it.
530    C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to
531    C disable this section of code.
532            IF (tempMultiDimAdvec) THEN
533              CALL GAD_ADVECTION(bi,bj,tempAdvScheme,GAD_TEMPERATURE,
534         U                      theta,gT,
535         I                      myTime,myIter,myThid)
536            ENDIF
537            IF (saltMultiDimAdvec) THEN
538              CALL GAD_ADVECTION(bi,bj,saltAdvScheme,GAD_SALINITY,
539         U                      salt,gS,
540         I                      myTime,myIter,myThid)
541            ENDIF
542    C Since passive tracers are configurable separately from T,S we
543    C call the multi-dimensional method for PTRACERS regardless
544    C of whether multiDimAdvection is set or not.
545    #ifdef ALLOW_PTRACERS
546            IF ( usePTRACERS ) THEN
547             CALL PTRACERS_ADVECTION( bi,bj,myIter,myTime,myThid )
548            ENDIF
549    #endif /* ALLOW_PTRACERS */
550    #endif /* DISABLE_MULTIDIM_ADVECTION */
551    
552  C--     Start of thermodynamics loop  C--     Start of thermodynamics loop
553          DO k=Nr,1,-1          DO k=Nr,1,-1
# Line 443  C--     Start of thermodynamics loop Line 555  C--     Start of thermodynamics loop
555  C? Patrick Is this formula correct?  C? Patrick Is this formula correct?
556  cph Yes, but I rewrote it.  cph Yes, but I rewrote it.
557  cph Also, the KappaR? need the index and subscript k!  cph Also, the KappaR? need the index and subscript k!
558           kkey = (ikey-1)*Nr + k           kkey = (itdkey-1)*Nr + k
559  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
560    
561  C--       km1    Points to level above k (=k-1)  C--       km1    Points to level above k (=k-1)
# Line 465  C--       Get temporary terms used by te Line 577  C--       Get temporary terms used by te
577       O         xA,yA,uTrans,vTrans,rTrans,maskUp,       O         xA,yA,uTrans,vTrans,rTrans,maskUp,
578       I         myThid)       I         myThid)
579    
580    #ifdef ALLOW_GMREDI
581    
582    C--   Residual transp = Bolus transp + Eulerian transp
583              IF (useGMRedi) THEN
584                CALL GMREDI_CALC_UVFLOW(
585         &            uTrans, vTrans, bi, bj, k, myThid)
586                IF (K.GE.2) CALL GMREDI_CALC_WFLOW(
587         &                    rTrans, bi, bj, k, myThid)
588              ENDIF
589    
590    #ifdef ALLOW_AUTODIFF_TAMC
591    #ifdef GM_BOLUS_ADVEC
592    CADJ STORE uTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
593    CADJ STORE vTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
594    CADJ STORE rTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
595    #endif
596    #endif /* ALLOW_AUTODIFF_TAMC */
597    
598    #endif /* ALLOW_GMREDI */
599    
600  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
601  CADJ STORE KappaRT(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE KappaRT(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte
602  CADJ STORE KappaRS(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE KappaRS(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte
# Line 492  C        and step forward storing result Line 624  C        and step forward storing result
624       I         xA,yA,uTrans,vTrans,rTrans,maskUp,       I         xA,yA,uTrans,vTrans,rTrans,maskUp,
625       I         KappaRT,       I         KappaRT,
626       U         fVerT,       U         fVerT,
627       I         myTime, myThid)       I         myTime,myIter,myThid)
            tauAB = 0.5d0 + abEps  
628             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
629       I         bi,bj,iMin,iMax,jMin,jMax,k,tauAB,       I         bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,
630       I         theta, gT,       I         theta, gT,
      U         gTnm1,  
631       I         myIter, myThid)       I         myIter, myThid)
632           ENDIF           ENDIF
633    cswdice ---- add ---
634    #ifdef ALLOW_THERM_SEAICE
635           if (k.eq.1) then
636            call ICE_FREEZE(bi, bj, iMin, iMax, jMin, jMax, myThid )
637           endif
638    #endif
639    cswdice -- end add ---
640           IF ( saltStepping ) THEN           IF ( saltStepping ) THEN
641             CALL CALC_GS(             CALL CALC_GS(
642       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
643       I         xA,yA,uTrans,vTrans,rTrans,maskUp,       I         xA,yA,uTrans,vTrans,rTrans,maskUp,
644       I         KappaRS,       I         KappaRS,
645       U         fVerS,       U         fVerS,
646       I         myTime, myThid)       I         myTime,myIter,myThid)
            tauAB = 0.5d0 + abEps  
647             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
648       I         bi,bj,iMin,iMax,jMin,jMax,k,tauAB,       I         bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,
649       I         salt, gS,       I         salt, gS,
      U         gSnm1,  
650       I         myIter, myThid)       I         myIter, myThid)
651           ENDIF           ENDIF
652  #ifdef ALLOW_PASSIVE_TRACER  #ifdef ALLOW_PASSIVE_TRACER
# Line 521  C        and step forward storing result Line 656  C        and step forward storing result
656       I         xA,yA,uTrans,vTrans,rTrans,maskUp,       I         xA,yA,uTrans,vTrans,rTrans,maskUp,
657       I         KappaRT,       I         KappaRT,
658       U         fVerTr1,       U         fVerTr1,
659       I         myTime, myThid)       I         myTime,myIter,myThid)
            tauAB = 0.5d0 + abEps  
660             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
661       I         bi,bj,iMin,iMax,jMin,jMax,k,tauAB,       I         bi,bj,iMin,iMax,jMin,jMax,k,tracerAdvScheme,
662       I         Tr1, gTr1,       I         Tr1, gTr1,
663       U         gTr1NM1,       I         myIter,myThid)
      I         myIter, myThid)  
664           ENDIF           ENDIF
665  #endif  #endif
666    #ifdef ALLOW_PTRACERS
667             IF ( usePTRACERS ) THEN
668               CALL PTRACERS_INTEGERATE(
669         I         bi,bj,k,
670         I         xA,yA,uTrans,vTrans,rTrans,maskUp,
671         X         KappaRS,
672         I         myIter,myTime,myThid)
673             ENDIF
674    #endif /* ALLOW_PTRACERS */
675    
676  #ifdef   ALLOW_OBCS  #ifdef   ALLOW_OBCS
677  C--      Apply open boundary conditions  C--      Apply open boundary conditions
678           IF (useOBCS) THEN           IF (useOBCS) THEN
679             CALL OBCS_APPLY_TS( bi, bj, k, gTnm1, gSnm1, myThid )             CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )
680           END IF           END IF
681  #endif   /* ALLOW_OBCS */  #endif   /* ALLOW_OBCS */
682    
683  C--      Freeze water  C--      Freeze water
684           IF (allowFreezing) THEN           IF ( allowFreezing .AND. .NOT. useSEAICE ) THEN
685  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
686  CADJ STORE gTNm1(:,:,k,bi,bj) = comlev1_bibj_k  CADJ STORE gT(:,:,k,bi,bj) = comlev1_bibj_k
687  CADJ &   , key = kkey, byte = isbyte  CADJ &   , key = kkey, byte = isbyte
688  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
689              CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )              CALL FREEZE( bi, bj, iMin, iMax, jMin, jMax, k, myThid )
# Line 550  CADJ &   , key = kkey, byte = isbyte Line 692  CADJ &   , key = kkey, byte = isbyte
692  C--     end of thermodynamic k loop (Nr:1)  C--     end of thermodynamic k loop (Nr:1)
693          ENDDO          ENDDO
694    
695    cswdice -- add ---
696    #ifdef ALLOW_THERM_SEAICE
697    c timeaveraging for ice model values
698               CALL ICE_AVE(bi,bj,iMin,iMax,jMin,jMax,myThid )
699    #endif
700    cswdice --- end add ---
701    
702    
703    
 #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 */  
704    
705  C--     Implicit diffusion  C--     Implicit diffusion
706          IF (implicitDiffusion) THEN          IF (implicitDiffusion) THEN
707    
708           IF (tempStepping) THEN           IF (tempStepping) THEN
709  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
710              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  
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, KappaRT, recip_HFacC,       I         deltaTtracer, KappaRT, recip_HFacC,
715       U         gTNm1,       U         gT,
716       I         myThid )       I         myThid )
717           ENDIF           ENDIF
718    
719           IF (saltStepping) THEN           IF (saltStepping) THEN
720  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
721           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  
722  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
723              CALL IMPLDIFF(              CALL IMPLDIFF(
724       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
725       I         deltaTtracer, KappaRS, recip_HFacC,       I         deltaTtracer, KappaRS, recip_HFacC,
726       U         gSNm1,       U         gS,
727       I         myThid )       I         myThid )
728           ENDIF           ENDIF
729    
730  #ifdef ALLOW_PASSIVE_TRACER  #ifdef ALLOW_PASSIVE_TRACER
731           IF (tr1Stepping) THEN           IF (tr1Stepping) THEN
732  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
733  CADJ STORE gTr1Nm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE gTr1(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
734  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
735            CALL IMPLDIFF(            CALL IMPLDIFF(
736       I      bi, bj, iMin, iMax, jMin, jMax,       I      bi, bj, iMin, iMax, jMin, jMax,
737       I      deltaTtracer, KappaRT, recip_HFacC,       I      deltaTtracer, KappaRT, recip_HFacC,
738       U      gTr1Nm1,       U      gTr1,
739       I      myThid )       I      myThid )
740           ENDIF           ENDIF
741  #endif  #endif
742    
743    #ifdef ALLOW_PTRACERS
744    C Vertical diffusion (implicit) for passive tracers
745             IF ( usePTRACERS ) THEN
746               CALL PTRACERS_IMPLDIFF( bi,bj,KappaRS,myThid )
747             ENDIF
748    #endif /* ALLOW_PTRACERS */
749    
750  #ifdef   ALLOW_OBCS  #ifdef   ALLOW_OBCS
751  C--      Apply open boundary conditions  C--      Apply open boundary conditions
752           IF (useOBCS) THEN           IF (useOBCS) THEN
753             DO K=1,Nr             DO K=1,Nr
754               CALL OBCS_APPLY_TS( bi, bj, k, gTnm1, gSnm1, myThid )               CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )
755             ENDDO             ENDDO
756           END IF           END IF
757  #endif   /* ALLOW_OBCS */  #endif   /* ALLOW_OBCS */
# Line 613  C--      Apply open boundary conditions Line 759  C--      Apply open boundary conditions
759  C--     End If implicitDiffusion  C--     End If implicitDiffusion
760          ENDIF          ENDIF
761    
762    #endif /* SINGLE_LAYER_MODE */
763    
764  Ccs-  Ccs-
765         ENDDO         ENDDO
766        ENDDO        ENDDO
767    
768  #ifdef ALLOW_AIM  #ifdef ALLOW_AIM
769        IF ( useAIM ) THEN  c     IF ( useAIM ) THEN
770         CALL AIM_AIM2DYN_EXCHANGES( myTime, myThid )  c      CALL AIM_AIM2DYN_EXCHANGES( myTime, myThid )
771        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  
772  #endif /* ALLOW_AIM */  #endif /* ALLOW_AIM */
773    c     IF ( staggerTimeStep ) THEN
774    c      IF ( useAIM .OR. useCubedSphereExchange ) THEN
775    c        IF (tempStepping) _EXCH_XYZ_R8(gT,myThid)
776    c        IF (saltStepping) _EXCH_XYZ_R8(gS,myThid)
777    c      ELSEIF ( useGMRedi .AND. Oly.LT.4 ) THEN
778    cc       .AND. GM_AdvForm .AND. .NOT.GM_AdvSeparate ) THEN
779    c        IF (tempMultiDimAdvec) _EXCH_XYZ_R8(gT,myThid)
780    c        IF (saltMultiDimAdvec) _EXCH_XYZ_R8(gS,myThid)
781    c      ENDIF  
782    c     ENDIF  
783    
784    #ifndef DISABLE_DEBUGMODE
785          If (debugMode) THEN
786           CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (THERMODYNAMICS)',myThid)
787           CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (THERMODYNAMICS)',myThid)
788           CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (THERMODYNAMICS)',myThid)
789           CALL DEBUG_STATS_RL(Nr,theta,'Theta (THERMODYNAMICS)',myThid)
790           CALL DEBUG_STATS_RL(Nr,salt,'Salt (THERMODYNAMICS)',myThid)
791           CALL DEBUG_STATS_RL(Nr,Gt,'Gt (THERMODYNAMICS)',myThid)
792           CALL DEBUG_STATS_RL(Nr,Gs,'Gs (THERMODYNAMICS)',myThid)
793           CALL DEBUG_STATS_RL(Nr,GtNm1,'GtNm1 (THERMODYNAMICS)',myThid)
794           CALL DEBUG_STATS_RL(Nr,GsNm1,'GsNm1 (THERMODYNAMICS)',myThid)
795    #ifdef ALLOW_PTRACERS
796           IF ( usePTRACERS ) THEN
797             CALL PTRACERS_DEBUG(myThid)
798           ENDIF
799    #endif /* ALLOW_PTRACERS */
800          ENDIF
801    #endif
802    
803        RETURN        RETURN
804        END        END
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.35
  ViewVC Help
Powered by ViewVC 1.1.22