/[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.37 by jmc, Fri Feb 7 21:59:53 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 59  C                                      i Line 123  C                                      i
123  C                                      so we need an fVer for each  C                                      so we need an fVer for each
124  C                                      variable.  C                                      variable.
125  C     rhoK, rhoKM1   - Density at current level, and level above  C     rhoK, rhoKM1   - Density at current level, and level above
 C     phiHyd         - Hydrostatic part of the potential phiHydi.  
 C                      In z coords phiHydiHyd is the hydrostatic  
 C                      Potential (=pressure/rho0) anomaly  
 C                      In p coords phiHydiHyd is the geopotential  
 C                      surface height anomaly.  
126  C     phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean)  C     phiSurfX, - gradient of Surface potentiel (Pressure/rho, ocean)
127  C     phiSurfY             or geopotentiel (atmos) in X and Y direction  C     phiSurfY             or geopotentiel (atmos) in X and Y direction
128  C     KappaRT,       - Total diffusion in vertical for T and S.  C     KappaRT,       - Total diffusion in vertical for T and S.
# Line 74  C     bi, bj Line 133  C     bi, bj
133  C     k, kup,        - Index for layer above and below. kup and kDown  C     k, kup,        - Index for layer above and below. kup and kDown
134  C     kDown, km1       are switched with layer to be the appropriate  C     kDown, km1       are switched with layer to be the appropriate
135  C                      index into fVerTerm.  C                      index into fVerTerm.
 C     tauAB - Adams-Bashforth timestepping weight: 0=forward ; 1/2=Adams-Bashf.  
136        _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
137        _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
138        _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 142  C     tauAB - Adams-Bashforth timesteppi
142        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerT   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
143        _RL fVerS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerS   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
144        _RL fVerTr1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerTr1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
       _RL phiHyd  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
145        _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhokm1  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
146        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rhok    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
147        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL phiSurfX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 94  C     tauAB - Adams-Bashforth timesteppi Line 151  C     tauAB - Adams-Bashforth timesteppi
151        _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaX  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
152        _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaY  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
153        _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL sigmaR  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
154        _RL tauAB  C     This is currently used by IVDC and Diagnostics
   
 C This is currently used by IVDC and Diagnostics  
155        _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL ConvectCount (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
   
156        INTEGER iMin, iMax        INTEGER iMin, iMax
157        INTEGER jMin, jMax        INTEGER jMin, jMax
158        INTEGER bi, bj        INTEGER bi, bj
159        INTEGER i, j        INTEGER i, j
160        INTEGER k, km1, kup, kDown        INTEGER k, km1, kup, kDown
161    
162  Cjmc : add for phiHyd output <- but not working if multi tile per CPU  CEOP
 c     CHARACTER*(MAX_LEN_MBUF) suff  
 c     LOGICAL  DIFFERENT_MULTIPLE  
 c     EXTERNAL DIFFERENT_MULTIPLE  
 Cjmc(end)  
163    
 C---    The algorithm...  
 C  
 C       "Correction Step"  
 C       =================  
 C       Here we update the horizontal velocities with the surface  
 C       pressure such that the resulting flow is either consistent  
 C       with the free-surface evolution or the rigid-lid:  
 C         U[n] = U* + dt x d/dx P  
 C         V[n] = V* + dt x d/dy P  
 C  
 C       "Calculation of Gs"  
 C       ===================  
 C       This is where all the accelerations and tendencies (ie.  
 C       physics, parameterizations etc...) are calculated  
 C         rho = rho ( theta[n], salt[n] )  
 C         b   = b(rho, theta)  
 C         K31 = K31 ( rho )  
 C         Gu[n] = Gu( u[n], v[n], wVel, b, ... )  
 C         Gv[n] = Gv( u[n], v[n], wVel, b, ... )  
 C         Gt[n] = Gt( theta[n], u[n], v[n], wVel, K31, ... )  
 C         Gs[n] = Gs( salt[n], u[n], v[n], wVel, K31, ... )  
 C  
 C       "Time-stepping" or "Prediction"  
 C       ================================  
 C       The models variables are stepped forward with the appropriate  
 C       time-stepping scheme (currently we use Adams-Bashforth II)  
 C       - For momentum, the result is always *only* a "prediction"  
 C       in that the flow may be divergent and will be "corrected"  
 C       later with a surface pressure gradient.  
 C       - Normally for tracers the result is the new field at time  
 C       level [n+1} *BUT* in the case of implicit diffusion the result  
 C       is also *only* a prediction.  
 C       - We denote "predictors" with an asterisk (*).  
 C         U* = U[n] + dt x ( 3/2 Gu[n] - 1/2 Gu[n-1] )  
 C         V* = V[n] + dt x ( 3/2 Gv[n] - 1/2 Gv[n-1] )  
 C         theta[n+1] = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )  
 C         salt[n+1] = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )  
 C       With implicit diffusion:  
 C         theta* = theta[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )  
 C         salt* = salt[n] + dt x ( 3/2 Gt[n] - 1/2 atG[n-1] )  
 C         (1 + dt * K * d_zz) theta[n] = theta*  
 C         (1 + dt * K * d_zz) salt[n] = salt*  
 C---  
   
164  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
165  C--   dummy statement to end declaration part  C--   dummy statement to end declaration part
166        ikey = 1        ikey = 1
167          itdkey = 1
168  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
169    
170  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 178  C     uninitialised but inert locations.
178          yA(i,j)      = 0. _d 0          yA(i,j)      = 0. _d 0
179          uTrans(i,j)  = 0. _d 0          uTrans(i,j)  = 0. _d 0
180          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  
181          rhok   (i,j) = 0. _d 0          rhok   (i,j) = 0. _d 0
182          phiSurfX(i,j) = 0. _d 0          phiSurfX(i,j) = 0. _d 0
183          phiSurfY(i,j) = 0. _d 0          phiSurfY(i,j) = 0. _d 0
# Line 195  CHPF$ INDEPENDENT Line 195  CHPF$ INDEPENDENT
195  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
196  C--    HPF directive to help TAMC  C--    HPF directive to help TAMC
197  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS  CHPF$  INDEPENDENT, NEW (rTrans,fVerT,fVerS
198  CHPF$&                  ,phiHyd,utrans,vtrans,xA,yA  CHPF$&                  ,utrans,vtrans,xA,yA
199  CHPF$&                  ,KappaRT,KappaRS  CHPF$&                  ,KappaRT,KappaRS
200  CHPF$&                  )  CHPF$&                  )
201  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
# Line 205  CHPF$&                  ) Line 205  CHPF$&                  )
205  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
206            act1 = bi - myBxLo(myThid)            act1 = bi - myBxLo(myThid)
207            max1 = myBxHi(myThid) - myBxLo(myThid) + 1            max1 = myBxHi(myThid) - myBxLo(myThid) + 1
   
208            act2 = bj - myByLo(myThid)            act2 = bj - myByLo(myThid)
209            max2 = myByHi(myThid) - myByLo(myThid) + 1            max2 = myByHi(myThid) - myByLo(myThid) + 1
   
210            act3 = myThid - 1            act3 = myThid - 1
211            max3 = nTx*nTy            max3 = nTx*nTy
   
212            act4 = ikey_dynamics - 1            act4 = ikey_dynamics - 1
213              itdkey = (act1 + 1) + act2*max1
           ikey = (act1 + 1) + act2*max1  
214       &                      + act3*max1*max2       &                      + act3*max1*max2
215       &                      + act4*max1*max2*max3       &                      + act4*max1*max2*max3
216  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
# Line 229  C--     Set up work arrays that need val Line 225  C--     Set up work arrays that need val
225            fVerS  (i,j,2) = 0. _d 0            fVerS  (i,j,2) = 0. _d 0
226            fVerTr1(i,j,1) = 0. _d 0            fVerTr1(i,j,1) = 0. _d 0
227            fVerTr1(i,j,2) = 0. _d 0            fVerTr1(i,j,2) = 0. _d 0
228              rhoKM1 (i,j)   = 0. _d 0
229           ENDDO           ENDDO
230          ENDDO          ENDDO
231    
# Line 236  C--     Set up work arrays that need val Line 233  C--     Set up work arrays that need val
233           DO j=1-OLy,sNy+OLy           DO j=1-OLy,sNy+OLy
234            DO i=1-OLx,sNx+OLx            DO i=1-OLx,sNx+OLx
235  C This is currently also used by IVDC and Diagnostics  C This is currently also used by IVDC and Diagnostics
236               sigmaX(i,j,k) = 0. _d 0
237               sigmaY(i,j,k) = 0. _d 0
238               sigmaR(i,j,k) = 0. _d 0
239             ConvectCount(i,j,k) = 0.             ConvectCount(i,j,k) = 0.
240             KappaRT(i,j,k) = 0. _d 0             KappaRT(i,j,k)    = 0. _d 0
241             KappaRS(i,j,k) = 0. _d 0             KappaRS(i,j,k)    = 0. _d 0
242    #ifdef ALLOW_AUTODIFF_TAMC
243    cph all the following init. are necessary for TAF
244    cph although some of these are re-initialised later.
245               gT(i,j,k,bi,bj)   = 0. _d 0
246               gS(i,j,k,bi,bj)   = 0. _d 0
247    # ifdef ALLOW_PASSIVE_TRACER
248               gTr1(i,j,k,bi,bj) = 0. _d 0
249    # endif
250    # ifdef ALLOW_GMREDI
251               Kwx(i,j,k,bi,bj)  = 0. _d 0
252               Kwy(i,j,k,bi,bj)  = 0. _d 0
253               Kwz(i,j,k,bi,bj)  = 0. _d 0
254    #  ifdef GM_NON_UNITY_DIAGONAL
255               Kux(i,j,k,bi,bj)  = 0. _d 0
256               Kvy(i,j,k,bi,bj)  = 0. _d 0
257    #  endif
258    #  ifdef GM_EXTRA_DIAGONAL
259               Kuz(i,j,k,bi,bj)  = 0. _d 0
260               Kvz(i,j,k,bi,bj)  = 0. _d 0
261    #  endif
262    #  ifdef GM_BOLUS_ADVEC
263               GM_PsiX(i,j,k,bi,bj)  = 0. _d 0
264               GM_PsiY(i,j,k,bi,bj)  = 0. _d 0
265    #  endif
266    #  ifdef GM_VISBECK_VARIABLE_K
267               VisbeckK(i,j,bi,bj)   = 0. _d 0
268    #  endif
269    # endif /* ALLOW_GMREDI */
270    #endif /* ALLOW_AUTODIFF_TAMC */
271            ENDDO            ENDDO
272           ENDDO           ENDDO
273          ENDDO          ENDDO
274    
275          iMin = 1-OLx+1          iMin = 1-OLx
276          iMax = sNx+OLx          iMax = sNx+OLx
277          jMin = 1-OLy+1          jMin = 1-OLy
278          jMax = sNy+OLy          jMax = sNy+OLy
279    
280    
281  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
282  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
283  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
284    #ifdef ALLOW_KPP
285    CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
286    CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
287    #endif
288  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
289    
290  C--     Start of diagnostic loop  C--     Start of diagnostic loop
# Line 262  C? Patrick, is this formula correct now Line 295  C? Patrick, is this formula correct now
295  C? Do we still need this?  C? Do we still need this?
296  cph kkey formula corrected.  cph kkey formula corrected.
297  cph Needed for rhok, rhokm1, in the case useGMREDI.  cph Needed for rhok, rhokm1, in the case useGMREDI.
298           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  
299  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
300    
301  C--       Integrate continuity vertically for vertical velocity  C--       Integrate continuity vertically for vertical velocity
302            CALL INTEGRATE_FOR_W(  c         CALL INTEGRATE_FOR_W(
303       I                         bi, bj, k, uVel, vVel,  c    I                         bi, bj, k, uVel, vVel,
304       O                         wVel,  c    O                         wVel,
305       I                         myThid )  c    I                         myThid )
306    
307  #ifdef    ALLOW_OBCS  #ifdef    ALLOW_OBCS
308  #ifdef    ALLOW_NONHYDROSTATIC  #ifdef    ALLOW_NONHYDROSTATIC
309  C--       Apply OBC to W if in N-H mode  C--       Apply OBC to W if in N-H mode
310            IF (useOBCS.AND.nonHydrostatic) THEN  c         IF (useOBCS.AND.nonHydrostatic) THEN
311              CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )  c           CALL OBCS_APPLY_W( bi, bj, k, wVel, myThid )
312            ENDIF  c         ENDIF
313  #endif    /* ALLOW_NONHYDROSTATIC */  #endif    /* ALLOW_NONHYDROSTATIC */
314  #endif    /* ALLOW_OBCS */  #endif    /* ALLOW_OBCS */
315    
316    C--       Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h
317    C--       MOST of THERMODYNAMICS will be disabled
318    #ifndef SINGLE_LAYER_MODE
319    
320  C--       Calculate gradients of potential density for isoneutral  C--       Calculate gradients of potential density for isoneutral
321  C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)  C         slope terms (e.g. GM/Redi tensor or IVDC diffusivity)
322  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 324  c         IF ( k.GT.1 .AND. (useGMRedi.O
324  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
325  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
326  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
327    CADJ STORE pressure(:,:,k,bi,bj) =
328    CADJ &     comlev1_bibj_k, key=kkey, byte=isbyte
329  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
330              CALL FIND_RHO(              CALL FIND_RHO(
331       I        bi, bj, iMin, iMax, jMin, jMax, k, k, eosType,       I        bi, bj, iMin, iMax, jMin, jMax, k, k,
332       I        theta, salt,       I        theta, salt,
333       O        rhoK,       O        rhoK,
334       I        myThid )       I        myThid )
335    
336              IF (k.GT.1) THEN              IF (k.GT.1) THEN
337  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
338  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
339  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
340    CADJ STORE pressure(:,:,k-1,bi,bj) =
341    CADJ &     comlev1_bibj_k, key=kkey, byte=isbyte
342  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
343               CALL FIND_RHO(               CALL FIND_RHO(
344       I        bi, bj, iMin, iMax, jMin, jMax, k-1, k, eosType,       I        bi, bj, iMin, iMax, jMin, jMax, k-1, k,
345       I        theta, salt,       I        theta, salt,
346       O        rhoKm1,       O        rhoKm1,
347       I        myThid )       I        myThid )
# Line 313  CADJ STORE salt (:,:,k-1,bi,bj) = comlev Line 353  CADJ STORE salt (:,:,k-1,bi,bj) = comlev
353       I             myThid )       I             myThid )
354            ENDIF            ENDIF
355    
356    #ifdef ALLOW_AUTODIFF_TAMC
357    CADJ STORE rhok   (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte
358    CADJ STORE rhokm1 (:,:) = comlev1_bibj_k ,       key=kkey, byte=isbyte
359    #endif /* ALLOW_AUTODIFF_TAMC */
360  C--       Implicit Vertical Diffusion for Convection  C--       Implicit Vertical Diffusion for Convection
361  c ==> should use sigmaR !!!  c ==> should use sigmaR !!!
362            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 367  c ==> should use sigmaR !!!
367       I        myTime, myIter, myThid)       I        myTime, myIter, myThid)
368            ENDIF            ENDIF
369    
370    #endif /* SINGLE_LAYER_MODE */
371    
372  C--     end of diagnostic k loop (Nr:1)  C--     end of diagnostic k loop (Nr:1)
373          ENDDO          ENDDO
374    
375  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
376  cph avoids recomputation of integrate_for_w  cph avoids recomputation of integrate_for_w
377  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE wvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
378    CADJ STORE pressure (:,:,:,bi,bj) =
379    CADJ &     comlev1_bibj, key=itdkey, byte=isbyte
380  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
381    
382  #ifdef  ALLOW_OBCS  #ifdef  ALLOW_OBCS
383  C--     Calculate future values on open boundaries  C--     Calculate future values on open boundaries
384          IF (useOBCS) THEN          IF (useOBCS) THEN
385            CALL OBCS_CALC( bi, bj, myTime+deltaT,            CALL OBCS_CALC( bi, bj, myTime+deltaT, myIter+1,
386       I            uVel, vVel, wVel, theta, salt,       I            uVel, vVel, wVel, theta, salt,
387       I            myThid )       I            myThid )
388          ENDIF          ENDIF
389  #endif  /* ALLOW_OBCS */  #endif  /* ALLOW_OBCS */
390    
391    
392    c********************************************
393    cswdice --- add ---
394    #ifdef ALLOW_THERM_SEAICE
395    C--     Determines forcing terms based on external fields
396    c--     including effects from ice
397            CALL ICE_FORCING(
398         I             bi, bj, iMin, iMax, jMin, jMax,
399         I             myThid )
400    #else
401    
402    cswdice --- end add ---
403    
404  C--     Determines forcing terms based on external fields  C--     Determines forcing terms based on external fields
405  C       relaxation terms, etc.  C       relaxation terms, etc.
406          CALL EXTERNAL_FORCING_SURF(          CALL EXTERNAL_FORCING_SURF(
407       I             bi, bj, iMin, iMax, jMin, jMax,       I             bi, bj, iMin, iMax, jMin, jMax,
408       I             myThid )       I             myThid )
409    cswdice --- add ----
410    #endif
411    cswdice --- end add ---
412    c******************************************
413    
414    
415    
416    
417  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
418  cph needed for KPP  cph needed for KPP
419  CADJ STORE surfacetendencyU(:,:,bi,bj)  CADJ STORE surfacetendencyU(:,:,bi,bj)
420  CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
421  CADJ STORE surfacetendencyV(:,:,bi,bj)  CADJ STORE surfacetendencyV(:,:,bi,bj)
422  CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
423  CADJ STORE surfacetendencyS(:,:,bi,bj)  CADJ STORE surfacetendencyS(:,:,bi,bj)
424  CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
425  CADJ STORE surfacetendencyT(:,:,bi,bj)  CADJ STORE surfacetendencyT(:,:,bi,bj)
426  CADJ &     = comlev1_bibj, key=ikey, byte=isbyte  CADJ &     = comlev1_bibj, key=itdkey, byte=isbyte
427  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
428    
429    C--     Attention: by defining "SINGLE_LAYER_MODE" in CPP_OPTIONS.h
430    C--     MOST of THERMODYNAMICS will be disabled
431    #ifndef SINGLE_LAYER_MODE
432    
433  #ifdef  ALLOW_GMREDI  #ifdef  ALLOW_GMREDI
434    
435  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
436  CADJ STORE sigmaX(:,:,:) = comlev1, key=ikey, byte=isbyte  cph storing here is needed only for one GMREDI_OPTIONS:
437  CADJ STORE sigmaY(:,:,:) = comlev1, key=ikey, byte=isbyte  cph define GM_BOLUS_ADVEC
438  CADJ STORE sigmaR(:,:,:) = comlev1, key=ikey, byte=isbyte  cph but I've avoided the #ifdef for now, in case more things change
439    CADJ STORE sigmaX(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte
440    CADJ STORE sigmaY(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte
441    CADJ STORE sigmaR(:,:,:)        = comlev1_bibj, key=itdkey, byte=isbyte
442  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
443    
444  C--     Calculate iso-neutral slopes for the GM/Redi parameterisation  C--     Calculate iso-neutral slopes for the GM/Redi parameterisation
445          IF (useGMRedi) THEN          IF (useGMRedi) THEN
446            DO k=1,Nr            CALL GMREDI_CALC_TENSOR(
447              CALL GMREDI_CALC_TENSOR(       I             bi, bj, iMin, iMax, jMin, jMax,
      I             bi, bj, iMin, iMax, jMin, jMax, k,  
448       I             sigmaX, sigmaY, sigmaR,       I             sigmaX, sigmaY, sigmaR,
449       I             myThid )       I             myThid )
           ENDDO  
450  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
451          ELSE          ELSE
452            DO k=1, Nr            CALL GMREDI_CALC_TENSOR_DUMMY(
453              CALL GMREDI_CALC_TENSOR_DUMMY(       I             bi, bj, iMin, iMax, jMin, jMax,
      I             bi, bj, iMin, iMax, jMin, jMax, k,  
454       I             sigmaX, sigmaY, sigmaR,       I             sigmaX, sigmaY, sigmaR,
455       I             myThid )       I             myThid )
           ENDDO  
456  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
457          ENDIF          ENDIF
458    
459  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
460  CADJ STORE Kwx(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE Kwx(:,:,:,bi,bj)   = comlev1_bibj, key=itdkey, byte=isbyte
461  CADJ STORE Kwy(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE Kwy(:,:,:,bi,bj)   = comlev1_bibj, key=itdkey, byte=isbyte
462  CADJ STORE Kwz(:,:,:,bi,bj)   = comlev1_bibj, key=ikey, byte=isbyte  CADJ STORE Kwz(:,:,:,bi,bj)   = comlev1_bibj, key=itdkey, byte=isbyte
463  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
464    
465  #endif  /* ALLOW_GMREDI */  #endif  /* ALLOW_GMREDI */
# Line 405  C--     Compute KPP mixing coefficients Line 478  C--     Compute KPP mixing coefficients
478    
479  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
480  CADJ STORE KPPghat   (:,:,:,bi,bj)  CADJ STORE KPPghat   (:,:,:,bi,bj)
 CADJ &   , KPPviscAz (:,:,:,bi,bj)  
481  CADJ &   , KPPdiffKzT(:,:,:,bi,bj)  CADJ &   , KPPdiffKzT(:,:,:,bi,bj)
482  CADJ &   , KPPdiffKzS(:,:,:,bi,bj)  CADJ &   , KPPdiffKzS(:,:,:,bi,bj)
483  CADJ &   , KPPfrac   (:,:  ,bi,bj)  CADJ &   , KPPfrac   (:,:  ,bi,bj)
484  CADJ &                 = comlev1_bibj, key=ikey, byte=isbyte  CADJ &                 = comlev1_bibj, key=itdkey, byte=isbyte
485  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
486    
487  #endif  /* ALLOW_KPP */  #endif  /* ALLOW_KPP */
488    
489  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
490  CADJ STORE KappaRT(:,:,:)     = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE KappaRT(:,:,:)     = comlev1_bibj, key=itdkey, byte=isbyte
491  CADJ STORE KappaRS(:,:,:)     = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE KappaRS(:,:,:)     = comlev1_bibj, key=itdkey, byte=isbyte
492  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE theta(:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
493  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE salt (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
494  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE uvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
495  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE vvel (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
496  #ifdef ALLOW_PASSIVE_TRACER  #ifdef ALLOW_PASSIVE_TRACER
497  CADJ STORE tr1  (:,:,:,bi,bj) = comlev1_bibj, key = ikey, byte = isbyte  CADJ STORE tr1  (:,:,:,bi,bj) = comlev1_bibj, key=itdkey, byte=isbyte
498  #endif  #endif
499  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
500    
501  #ifdef ALLOW_AIM  #ifdef ALLOW_AIM
502  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  
503          IF ( useAIM ) THEN          IF ( useAIM ) THEN
504           CALL TIMER_START('AIM_DO_ATMOS_PHYS      [DYNAMICS]', myThid)           CALL TIMER_START('AIM_DO_PHYSICS   [THERMODYNAMICS]', myThid)
505           CALL AIM_DO_ATMOS_PHYSICS( phiHyd, bi, bj, myTime, myThid )           CALL AIM_DO_PHYSICS( bi, bj, myTime, myIter, myThid )
506           CALL TIMER_STOP ('AIM_DO_ATMOS_PHYS      [DYNAMICS]', myThid)           CALL TIMER_STOP( 'AIM_DO_PHYSICS   [THERMODYNAMICS]', myThid)
507          ENDIF          ENDIF
508  #endif /* ALLOW_AIM */  #endif /* ALLOW_AIM */
509    
510    #ifdef ALLOW_TIMEAVE
511            IF (taveFreq.GT.0. .AND. ivdc_kappa.NE.0.) THEN
512              CALL TIMEAVE_CUMULATE(ConvectCountTave, ConvectCount, Nr,
513         I                               deltaTclock, bi, bj, myThid)
514            ENDIF
515    #endif /* ALLOW_TIMEAVE */
516    
517    #ifndef DISABLE_MULTIDIM_ADVECTION
518    C--     Some advection schemes are better calculated using a multi-dimensional
519    C       method in the absence of any other terms and, if used, is done here.
520    C
521    C The CPP flag DISABLE_MULTIDIM_ADVECTION is currently unset in GAD_OPTIONS.h
522    C The default is to use multi-dimensinal advection for non-linear advection
523    C schemes. However, for the sake of efficiency of the adjoint it is necessary
524    C to be able to exclude this scheme to avoid excessive storage and
525    C recomputation. It *is* differentiable, if you need it.
526    C Edit GAD_OPTIONS.h and #define DISABLE_MULTIDIM_ADVECTION to
527    C disable this section of code.
528            IF (tempMultiDimAdvec) THEN
529              CALL GAD_ADVECTION(bi,bj,tempAdvScheme,GAD_TEMPERATURE,
530         U                      theta,gT,
531         I                      myTime,myIter,myThid)
532            ENDIF
533            IF (saltMultiDimAdvec) THEN
534              CALL GAD_ADVECTION(bi,bj,saltAdvScheme,GAD_SALINITY,
535         U                      salt,gS,
536         I                      myTime,myIter,myThid)
537            ENDIF
538    C Since passive tracers are configurable separately from T,S we
539    C call the multi-dimensional method for PTRACERS regardless
540    C of whether multiDimAdvection is set or not.
541    #ifdef ALLOW_PTRACERS
542            IF ( usePTRACERS ) THEN
543             CALL PTRACERS_ADVECTION( bi,bj,myIter,myTime,myThid )
544            ENDIF
545    #endif /* ALLOW_PTRACERS */
546    #endif /* DISABLE_MULTIDIM_ADVECTION */
547    
548  C--     Start of thermodynamics loop  C--     Start of thermodynamics loop
549          DO k=Nr,1,-1          DO k=Nr,1,-1
# Line 443  C--     Start of thermodynamics loop Line 551  C--     Start of thermodynamics loop
551  C? Patrick Is this formula correct?  C? Patrick Is this formula correct?
552  cph Yes, but I rewrote it.  cph Yes, but I rewrote it.
553  cph Also, the KappaR? need the index and subscript k!  cph Also, the KappaR? need the index and subscript k!
554           kkey = (ikey-1)*Nr + k           kkey = (itdkey-1)*Nr + k
555  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
556    
557  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 573  C--       Get temporary terms used by te
573       O         xA,yA,uTrans,vTrans,rTrans,maskUp,       O         xA,yA,uTrans,vTrans,rTrans,maskUp,
574       I         myThid)       I         myThid)
575    
576    #ifdef ALLOW_GMREDI
577    
578    C--   Residual transp = Bolus transp + Eulerian transp
579              IF (useGMRedi) THEN
580                CALL GMREDI_CALC_UVFLOW(
581         &            uTrans, vTrans, bi, bj, k, myThid)
582                IF (K.GE.2) CALL GMREDI_CALC_WFLOW(
583         &                    rTrans, bi, bj, k, myThid)
584              ENDIF
585    
586    #ifdef ALLOW_AUTODIFF_TAMC
587    #ifdef GM_BOLUS_ADVEC
588    CADJ STORE uTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
589    CADJ STORE vTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
590    CADJ STORE rTrans(:,:)    = comlev1_bibj_k, key=kkey, byte=isbyte
591    #endif
592    #endif /* ALLOW_AUTODIFF_TAMC */
593    
594    #endif /* ALLOW_GMREDI */
595    
596  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
597  CADJ STORE KappaRT(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte  CADJ STORE KappaRT(:,:,k)    = comlev1_bibj_k, key=kkey, byte=isbyte
598  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 620  C        and step forward storing result
620       I         xA,yA,uTrans,vTrans,rTrans,maskUp,       I         xA,yA,uTrans,vTrans,rTrans,maskUp,
621       I         KappaRT,       I         KappaRT,
622       U         fVerT,       U         fVerT,
623       I         myTime, myThid)       I         myTime,myIter,myThid)
            tauAB = 0.5d0 + abEps  
624             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
625       I         bi,bj,iMin,iMax,jMin,jMax,k,tauAB,       I         bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,
626       I         theta, gT,       I         theta, gT,
      U         gTnm1,  
627       I         myIter, myThid)       I         myIter, myThid)
628           ENDIF           ENDIF
629    cswdice ---- add ---
630    #ifdef ALLOW_THERM_SEAICE
631           if (k.eq.1) then
632            call ICE_FREEZE(bi, bj, iMin, iMax, jMin, jMax, myThid )
633           endif
634    #endif
635    cswdice -- end add ---
636           IF ( saltStepping ) THEN           IF ( saltStepping ) THEN
637             CALL CALC_GS(             CALL CALC_GS(
638       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,       I         bi,bj,iMin,iMax,jMin,jMax, k,km1,kup,kDown,
639       I         xA,yA,uTrans,vTrans,rTrans,maskUp,       I         xA,yA,uTrans,vTrans,rTrans,maskUp,
640       I         KappaRS,       I         KappaRS,
641       U         fVerS,       U         fVerS,
642       I         myTime, myThid)       I         myTime,myIter,myThid)
            tauAB = 0.5d0 + abEps  
643             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
644       I         bi,bj,iMin,iMax,jMin,jMax,k,tauAB,       I         bi,bj,iMin,iMax,jMin,jMax,k,saltAdvScheme,
645       I         salt, gS,       I         salt, gS,
      U         gSnm1,  
646       I         myIter, myThid)       I         myIter, myThid)
647           ENDIF           ENDIF
648  #ifdef ALLOW_PASSIVE_TRACER  #ifdef ALLOW_PASSIVE_TRACER
# Line 521  C        and step forward storing result Line 652  C        and step forward storing result
652       I         xA,yA,uTrans,vTrans,rTrans,maskUp,       I         xA,yA,uTrans,vTrans,rTrans,maskUp,
653       I         KappaRT,       I         KappaRT,
654       U         fVerTr1,       U         fVerTr1,
655       I         myTime, myThid)       I         myTime,myIter,myThid)
            tauAB = 0.5d0 + abEps  
656             CALL TIMESTEP_TRACER(             CALL TIMESTEP_TRACER(
657       I         bi,bj,iMin,iMax,jMin,jMax,k,tauAB,       I         bi,bj,iMin,iMax,jMin,jMax,k,tracerAdvScheme,
658       I         Tr1, gTr1,       I         Tr1, gTr1,
659       U         gTr1NM1,       I         myIter,myThid)
      I         myIter, myThid)  
660           ENDIF           ENDIF
661  #endif  #endif
662    #ifdef ALLOW_PTRACERS
663             IF ( usePTRACERS ) THEN
664               CALL PTRACERS_INTEGERATE(
665         I         bi,bj,k,
666         I         xA,yA,uTrans,vTrans,rTrans,maskUp,
667         X         KappaRS,
668         I         myIter,myTime,myThid)
669             ENDIF
670    #endif /* ALLOW_PTRACERS */
671    
672  #ifdef   ALLOW_OBCS  #ifdef   ALLOW_OBCS
673  C--      Apply open boundary conditions  C--      Apply open boundary conditions
674           IF (useOBCS) THEN           IF (useOBCS) THEN
675             CALL OBCS_APPLY_TS( bi, bj, k, gTnm1, gSnm1, myThid )             CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )
676           END IF           END IF
677  #endif   /* ALLOW_OBCS */  #endif   /* ALLOW_OBCS */
678    
679  C--      Freeze water  C--      Freeze water
680           IF (allowFreezing) THEN           IF ( allowFreezing .AND. .NOT. useSEAICE ) THEN
681  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
682  CADJ STORE gTNm1(:,:,k,bi,bj) = comlev1_bibj_k  CADJ STORE gT(:,:,k,bi,bj) = comlev1_bibj_k
683  CADJ &   , key = kkey, byte = isbyte  CADJ &   , key = kkey, byte = isbyte
684  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
685              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 688  CADJ &   , key = kkey, byte = isbyte
688  C--     end of thermodynamic k loop (Nr:1)  C--     end of thermodynamic k loop (Nr:1)
689          ENDDO          ENDDO
690    
691    cswdice -- add ---
692    #ifdef ALLOW_THERM_SEAICE
693    c timeaveraging for ice model values
694               CALL ICE_AVE(bi,bj,iMin,iMax,jMin,jMax,myThid )
695    #endif
696    cswdice --- end add ---
697    
698    
699    
 #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 */  
700    
701  C--     Implicit diffusion  C--     Implicit diffusion
702          IF (implicitDiffusion) THEN          IF (implicitDiffusion) THEN
703    
704           IF (tempStepping) THEN           IF (tempStepping) THEN
705  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
706              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  
707  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
708              CALL IMPLDIFF(              CALL IMPLDIFF(
709       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
710       I         deltaTtracer, KappaRT, recip_HFacC,       I         deltaTtracer, KappaRT, recip_HFacC,
711       U         gTNm1,       U         gT,
712       I         myThid )       I         myThid )
713           ENDIF           ENDIF
714    
715           IF (saltStepping) THEN           IF (saltStepping) THEN
716  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
717           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  
718  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
719              CALL IMPLDIFF(              CALL IMPLDIFF(
720       I         bi, bj, iMin, iMax, jMin, jMax,       I         bi, bj, iMin, iMax, jMin, jMax,
721       I         deltaTtracer, KappaRS, recip_HFacC,       I         deltaTtracer, KappaRS, recip_HFacC,
722       U         gSNm1,       U         gS,
723       I         myThid )       I         myThid )
724           ENDIF           ENDIF
725    
726  #ifdef ALLOW_PASSIVE_TRACER  #ifdef ALLOW_PASSIVE_TRACER
727           IF (tr1Stepping) THEN           IF (tr1Stepping) THEN
728  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
729  CADJ STORE gTr1Nm1(:,:,:,bi,bj) = comlev1_bibj , key=ikey, byte=isbyte  CADJ STORE gTr1(:,:,:,bi,bj) = comlev1_bibj , key=itdkey, byte=isbyte
730  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
731            CALL IMPLDIFF(            CALL IMPLDIFF(
732       I      bi, bj, iMin, iMax, jMin, jMax,       I      bi, bj, iMin, iMax, jMin, jMax,
733       I      deltaTtracer, KappaRT, recip_HFacC,       I      deltaTtracer, KappaRT, recip_HFacC,
734       U      gTr1Nm1,       U      gTr1,
735       I      myThid )       I      myThid )
736           ENDIF           ENDIF
737  #endif  #endif
738    
739    #ifdef ALLOW_PTRACERS
740    C Vertical diffusion (implicit) for passive tracers
741             IF ( usePTRACERS ) THEN
742               CALL PTRACERS_IMPLDIFF( bi,bj,KappaRS,myThid )
743             ENDIF
744    #endif /* ALLOW_PTRACERS */
745    
746  #ifdef   ALLOW_OBCS  #ifdef   ALLOW_OBCS
747  C--      Apply open boundary conditions  C--      Apply open boundary conditions
748           IF (useOBCS) THEN           IF (useOBCS) THEN
749             DO K=1,Nr             DO K=1,Nr
750               CALL OBCS_APPLY_TS( bi, bj, k, gTnm1, gSnm1, myThid )               CALL OBCS_APPLY_TS( bi, bj, k, gT, gS, myThid )
751             ENDDO             ENDDO
752           END IF           END IF
753  #endif   /* ALLOW_OBCS */  #endif   /* ALLOW_OBCS */
# Line 613  C--      Apply open boundary conditions Line 755  C--      Apply open boundary conditions
755  C--     End If implicitDiffusion  C--     End If implicitDiffusion
756          ENDIF          ENDIF
757    
758    #endif /* SINGLE_LAYER_MODE */
759    
760  Ccs-  Ccs-
761         ENDDO         ENDDO
762        ENDDO        ENDDO
763    
764  #ifdef ALLOW_AIM  #ifdef ALLOW_AIM
765        IF ( useAIM ) THEN  c     IF ( useAIM ) THEN
766         CALL AIM_AIM2DYN_EXCHANGES( myTime, myThid )  c      CALL AIM_AIM2DYN_EXCHANGES( myTime, myThid )
767        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  
768  #endif /* ALLOW_AIM */  #endif /* ALLOW_AIM */
769    c     IF ( staggerTimeStep ) THEN
770    c      IF ( useAIM .OR. useCubedSphereExchange ) THEN
771    c        IF (tempStepping) _EXCH_XYZ_R8(gT,myThid)
772    c        IF (saltStepping) _EXCH_XYZ_R8(gS,myThid)
773    c      ELSEIF ( useGMRedi .AND. Oly.LT.4 ) THEN
774    cc       .AND. GM_AdvForm .AND. .NOT.GM_AdvSeparate ) THEN
775    c        IF (tempMultiDimAdvec) _EXCH_XYZ_R8(gT,myThid)
776    c        IF (saltMultiDimAdvec) _EXCH_XYZ_R8(gS,myThid)
777    c      ENDIF  
778    c     ENDIF  
779    
780    #ifndef DISABLE_DEBUGMODE
781          If (debugMode) THEN
782           CALL DEBUG_STATS_RL(Nr,uVel,'Uvel (THERMODYNAMICS)',myThid)
783           CALL DEBUG_STATS_RL(Nr,vVel,'Vvel (THERMODYNAMICS)',myThid)
784           CALL DEBUG_STATS_RL(Nr,wVel,'Wvel (THERMODYNAMICS)',myThid)
785           CALL DEBUG_STATS_RL(Nr,theta,'Theta (THERMODYNAMICS)',myThid)
786           CALL DEBUG_STATS_RL(Nr,salt,'Salt (THERMODYNAMICS)',myThid)
787           CALL DEBUG_STATS_RL(Nr,Gt,'Gt (THERMODYNAMICS)',myThid)
788           CALL DEBUG_STATS_RL(Nr,Gs,'Gs (THERMODYNAMICS)',myThid)
789           CALL DEBUG_STATS_RL(Nr,GtNm1,'GtNm1 (THERMODYNAMICS)',myThid)
790           CALL DEBUG_STATS_RL(Nr,GsNm1,'GsNm1 (THERMODYNAMICS)',myThid)
791    #ifdef ALLOW_PTRACERS
792           IF ( usePTRACERS ) THEN
793             CALL PTRACERS_DEBUG(myThid)
794           ENDIF
795    #endif /* ALLOW_PTRACERS */
796          ENDIF
797    #endif
798    
799        RETURN        RETURN
800        END        END
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.37
  ViewVC Help
Powered by ViewVC 1.1.22