/[MITgcm]/MITgcm/pkg/generic_advdiff/gad_calc_rhs.F
ViewVC logotype

Diff of /MITgcm/pkg/generic_advdiff/gad_calc_rhs.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.25 by jmc, Fri Jun 25 18:19:20 2004 UTC revision 1.52 by jahn, Wed Apr 23 18:32:20 2008 UTC
# Line 7  CBOP Line 7  CBOP
7  C !ROUTINE: GAD_CALC_RHS  C !ROUTINE: GAD_CALC_RHS
8    
9  C !INTERFACE: ==========================================================  C !INTERFACE: ==========================================================
10        SUBROUTINE GAD_CALC_RHS(        SUBROUTINE GAD_CALC_RHS(
11       I           bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,       I           bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,
12       I           xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,       I           xA, yA, maskUp, uFld, vFld, wFld,
13       I           uVel, vVel, wVel,       I           uTrans, vTrans, rTrans, rTransKp1,
14       I           diffKh, diffK4, KappaRT, Tracer,       I           diffKh, diffK4, KappaR, TracerN, TracAB,
15       I           tracerIdentity, advectionScheme,       I           tracerIdentity, advectionScheme, vertAdvecScheme,
16       I           calcAdvection, implicitAdvection,       I           calcAdvection, implicitAdvection, applyAB_onTracer,
17         I           trUseGMRedi, trUseKPP,
18       U           fVerT, gTracer,       U           fVerT, gTracer,
19       I           myThid )       I           myTime, myIter, myThid )
20    
21  C !DESCRIPTION:  C !DESCRIPTION:
22  C Calculates the tendancy of a tracer due to advection and diffusion.  C Calculates the tendency of a tracer due to advection and diffusion.
23  C It calculates the fluxes in each direction indepentently and then  C It calculates the fluxes in each direction indepentently and then
24  C sets the tendancy to the divergence of these fluxes. The advective  C sets the tendency to the divergence of these fluxes. The advective
25  C fluxes are only calculated here when using the linear advection schemes  C fluxes are only calculated here when using the linear advection schemes
26  C otherwise only the diffusive and parameterized fluxes are calculated.  C otherwise only the diffusive and parameterized fluxes are calculated.
27  C  C
# Line 29  C \begin{equation*} Line 30  C \begin{equation*}
30  C {\bf F} = {\bf F}_{adv} + {\bf F}_{diff} +{\bf F}_{GM} + {\bf F}_{KPP}  C {\bf F} = {\bf F}_{adv} + {\bf F}_{diff} +{\bf F}_{GM} + {\bf F}_{KPP}
31  C \end{equation*}  C \end{equation*}
32  C  C
33  C The tendancy is the divergence of the fluxes:  C The tendency is the divergence of the fluxes:
34  C \begin{equation*}  C \begin{equation*}
35  C G_\theta = G_\theta + \nabla \cdot {\bf F}  C G_\theta = G_\theta + \nabla \cdot {\bf F}
36  C \end{equation*}  C \end{equation*}
37  C  C
38  C The tendancy is assumed to contain data on entry.  C The tendency is assumed to contain data on entry.
39    
40  C !USES: ===============================================================  C !USES: ===============================================================
41        IMPLICIT NONE        IMPLICIT NONE
# Line 54  C !INPUT PARAMETERS: =================== Line 55  C !INPUT PARAMETERS: ===================
55  C bi,bj            :: tile indices  C bi,bj            :: tile indices
56  C iMin,iMax        :: loop range for called routines  C iMin,iMax        :: loop range for called routines
57  C jMin,jMax        :: loop range for called routines  C jMin,jMax        :: loop range for called routines
58  C kup              :: index into 2 1/2D array, toggles between 1|2  C k                :: vertical index
59  C kdown            :: index into 2 1/2D array, toggles between 2|1  C kM1              :: =k-1 for k>1, =1 for k=1
60  C kp1              :: =k+1 for k<Nr, =Nr for k=Nr  C kUp              :: index into 2 1/2D array, toggles between 1|2
61    C kDown            :: index into 2 1/2D array, toggles between 2|1
62  C xA,yA            :: areas of X and Y face of tracer cells  C xA,yA            :: areas of X and Y face of tracer cells
63    C maskUp           :: 2-D array for mask at W points
64    C uFld,vFld,wFld   :: Local copy of velocity field (3 components)
65  C uTrans,vTrans    :: 2-D arrays of volume transports at U,V points  C uTrans,vTrans    :: 2-D arrays of volume transports at U,V points
66  C rTrans           :: 2-D arrays of volume transports at W points  C rTrans           :: 2-D arrays of volume transports at W points
67  C rTransKp1        :: 2-D array of volume trans at W pts, interf k+1  C rTransKp1        :: 2-D array of volume trans at W pts, interf k+1
 C maskUp           :: 2-D array for mask at W points  
 C uVel,vVel,wVel   :: 3 components of the velcity field (3-D array)  
68  C diffKh           :: horizontal diffusion coefficient  C diffKh           :: horizontal diffusion coefficient
69  C diffK4           :: bi-harmonic diffusion coefficient  C diffK4           :: bi-harmonic diffusion coefficient
70  C KappaRT          :: 3-D array for vertical diffusion coefficient  C KappaR           :: 2-D array for vertical diffusion coefficient, interf k
71  C Tracer           :: tracer field  C TracerN          :: tracer field @ time-step n (Note: only used
72    C                     if applying AB on tracer field rather than on tendency gTr)
73    C TracAB           :: current tracer field (@ time-step n if applying AB on gTr
74    C                     or extrapolated fwd in time to n+1/2 if applying AB on Tr)
75  C tracerIdentity   :: tracer identifier (required for KPP,GM)  C tracerIdentity   :: tracer identifier (required for KPP,GM)
76  C advectionScheme  :: advection scheme to use  C advectionScheme  :: advection scheme to use (Horizontal plane)
77    C vertAdvecScheme  :: advection scheme to use (Vertical direction)
78  C calcAdvection    :: =False if Advec computed with multiDim scheme  C calcAdvection    :: =False if Advec computed with multiDim scheme
79  C implicitAdvection:: =True if vertical Advec computed implicitly  C implicitAdvection:: =True if vertical Advec computed implicitly
80    C applyAB_onTracer :: apply Adams-Bashforth on Tracer (rather than on gTr)
81    C trUseGMRedi      :: true if this tracer uses GM-Redi
82    C trUseKPP         :: true if this tracer uses KPP
83    C myTime           :: current time
84    C myIter           :: iteration number
85  C myThid           :: thread number  C myThid           :: thread number
86        INTEGER bi,bj,iMin,iMax,jMin,jMax        INTEGER bi,bj,iMin,iMax,jMin,jMax
87        INTEGER k,kUp,kDown,kM1        INTEGER k,kUp,kDown,kM1
88        _RS xA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS xA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
89        _RS yA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS yA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
90          _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
91          _RL uFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
92          _RL vFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
93          _RL wFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
94        _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
95        _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
96        _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
97        _RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
       _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL uVel  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)  
       _RL vVel  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)  
       _RL wVel  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)  
98        _RL diffKh, diffK4        _RL diffKh, diffK4
99        _RL KappaRT(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL KappaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
100        _RL Tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL TracerN(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
101          _RL TracAB (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
102        INTEGER tracerIdentity        INTEGER tracerIdentity
103        INTEGER advectionScheme        INTEGER advectionScheme, vertAdvecScheme
104        LOGICAL calcAdvection        LOGICAL calcAdvection
105        LOGICAL implicitAdvection        LOGICAL implicitAdvection, applyAB_onTracer
106        INTEGER myThid        LOGICAL trUseGMRedi, trUseKPP
107          _RL     myTime
108          INTEGER myIter, myThid
109    
110  C !OUTPUT PARAMETERS: ==================================================  C !OUTPUT PARAMETERS: ==================================================
111  C gTracer          :: tendancy array  C gTracer          :: tendency array
112  C fVerT            :: 2 1/2D arrays for vertical advective flux  C fVerT            :: 2 1/2D arrays for vertical advective flux
113        _RL gTracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL gTracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
114        _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
# Line 103  C !LOCAL VARIABLES: ==================== Line 117  C !LOCAL VARIABLES: ====================
117  C i,j              :: loop indices  C i,j              :: loop indices
118  C df4              :: used for storing del^2 T for bi-harmonic term  C df4              :: used for storing del^2 T for bi-harmonic term
119  C fZon             :: zonal flux  C fZon             :: zonal flux
120  C fmer             :: meridional flux  C fMer             :: meridional flux
121  C af               :: advective flux  C af               :: advective flux
122  C df               :: diffusive flux  C df               :: diffusive flux
123  C localT           :: local copy of tracer field  C localT           :: local copy of tracer field
124    C locABT           :: local copy of (AB-extrapolated) tracer field
125    #ifdef ALLOW_DIAGNOSTICS
126          CHARACTER*8 diagName
127          CHARACTER*4 GAD_DIAG_SUFX, diagSufx
128          EXTERNAL    GAD_DIAG_SUFX
129    #endif
130        INTEGER i,j        INTEGER i,j
131        _RL df4   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL df4   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
132        _RL fZon  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fZon  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 114  C localT           :: local copy of trac Line 134  C localT           :: local copy of trac
134        _RL af    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL af    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
135        _RL df    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL df    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
136        _RL localT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL localT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
137          _RL locABT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
138        _RL advFac, rAdvFac        _RL advFac, rAdvFac
139    #ifdef GAD_SMOLARKIEWICZ_HACK
140          _RL outFlux, trac, fac, gTrFac
141    #endif
142  CEOP  CEOP
143    
144  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 123  C--   the kDown is still required Line 147  C--   the kDown is still required
147        fVerT(1,1,kDown) = fVerT(1,1,kDown)        fVerT(1,1,kDown) = fVerT(1,1,kDown)
148  #endif  #endif
149    
150    #ifdef ALLOW_DIAGNOSTICS
151    C--   Set diagnostic suffix for the current tracer
152          IF ( useDiagnostics ) THEN
153            diagSufx = GAD_DIAG_SUFX( tracerIdentity, myThid )
154          ENDIF
155    #endif
156    
157        advFac  = 0. _d 0        advFac  = 0. _d 0
158        IF (calcAdvection) advFac = 1. _d 0        IF (calcAdvection) advFac = 1. _d 0
159        rAdvFac = rkFac*advFac        rAdvFac = rkSign*advFac
160        IF (implicitAdvection) rAdvFac = 0. _d 0        IF (implicitAdvection) rAdvFac = 0. _d 0
161    
162        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
# Line 139  C--   the kDown is still required Line 170  C--   the kDown is still required
170        ENDDO        ENDDO
171    
172  C--   Make local copy of tracer array  C--   Make local copy of tracer array
173        DO j=1-OLy,sNy+OLy        IF ( applyAB_onTracer ) THEN
174         DO i=1-OLx,sNx+OLx          DO j=1-OLy,sNy+OLy
175          localT(i,j)=tracer(i,j,k,bi,bj)           DO i=1-OLx,sNx+OLx
176         ENDDO            localT(i,j)=TracerN(i,j,k,bi,bj)
177        ENDDO            locABT(i,j)= TracAB(i,j,k,bi,bj)
178             ENDDO
179            ENDDO
180          ELSE
181            DO j=1-OLy,sNy+OLy
182             DO i=1-OLx,sNx+OLx
183              localT(i,j)= TracAB(i,j,k,bi,bj)
184              locABT(i,j)= TracAB(i,j,k,bi,bj)
185             ENDDO
186            ENDDO
187          ENDIF
188    
189  C--   Unless we have already calculated the advection terms we initialize  C--   Unless we have already calculated the advection terms we initialize
190  C     the tendency to zero.  C     the tendency to zero.
# Line 172  C--   Initialize net flux in X direction Line 213  C--   Initialize net flux in X direction
213    
214  C-    Advective flux in X  C-    Advective flux in X
215        IF (calcAdvection) THEN        IF (calcAdvection) THEN
216        IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN          IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
217         CALL GAD_C2_ADV_X(bi,bj,k,uTrans,localT,af,myThid)            CALL GAD_C2_ADV_X(bi,bj,k,uTrans,locABT,af,myThid)
218        ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST
219         CALL GAD_FLUXLIMIT_ADV_X(       &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
220       &      bi,bj,k,deltaTtracer,uTrans,uVel,localT,af,myThid)            CALL GAD_DST2U1_ADV_X( bi,bj,k, advectionScheme, .TRUE.,
221        ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN       I            dTtracerLev(k), uTrans, uFld, locABT,
222         CALL GAD_U3_ADV_X(bi,bj,k,uTrans,localT,af,myThid)       O            af, myThid )
223        ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
224         CALL GAD_C4_ADV_X(bi,bj,k,uTrans,localT,af,myThid)            CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, .TRUE., dTtracerLev(k),
225        ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN       I            uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,
226         CALL GAD_DST3_ADV_X(       O            af, myThid )
227       &       bi,bj,k,deltaTtracer,uTrans,uVel,localT,af,myThid)          ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
228        ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN            CALL GAD_U3_ADV_X(bi,bj,k,uTrans,locABT,af,myThid)
229         CALL GAD_DST3FL_ADV_X(          ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
230       &       bi,bj,k,deltaTtracer,uTrans,uVel,localT,af,myThid)            CALL GAD_C4_ADV_X(bi,bj,k,uTrans,locABT,af,myThid)
231        ELSE          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
232         STOP 'GAD_CALC_RHS: Bad advectionScheme (X)'            CALL GAD_DST3_ADV_X( bi,bj,k, .TRUE., dTtracerLev(k),
233        ENDIF       I            uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,
234        DO j=1-Oly,sNy+Oly       O            af, myThid )
235         DO i=1-Olx,sNx+Olx          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
236          fZon(i,j) = fZon(i,j) + af(i,j)           IF ( inAdMode ) THEN
237         ENDDO  cph This block is to trick the adjoint:
238        ENDDO  cph IF inAdExact=.FALSE., we want to use DST3
239    cph with limiters in forward, but without limiters in reverse.
240              CALL GAD_DST3_ADV_X( bi,bj,k, .TRUE., dTtracerLev(k),
241         I           uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,
242         O           af, myThid )
243             ELSE
244              CALL GAD_DST3FL_ADV_X( bi,bj,k, .TRUE., dTtracerLev(k),
245         I           uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,
246         O           af, myThid )
247             ENDIF
248            ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN
249              CALL GAD_OS7MP_ADV_X( bi,bj,k, .TRUE., dTtracerLev(k),
250         I            uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,
251         O            af, myThid )
252            ELSE
253             STOP 'GAD_CALC_RHS: Bad advectionScheme (X)'
254            ENDIF
255            DO j=1-Oly,sNy+Oly
256             DO i=1-Olx,sNx+Olx
257              fZon(i,j) = fZon(i,j) + af(i,j)
258             ENDDO
259            ENDDO
260    #ifdef ALLOW_DIAGNOSTICS
261            IF ( useDiagnostics ) THEN
262              diagName = 'ADVx'//diagSufx
263              CALL DIAGNOSTICS_FILL(af,diagName, k,1, 2,bi,bj, myThid)
264            ENDIF
265    #endif
266        ENDIF        ENDIF
267    
268  C-    Diffusive flux in X  C-    Diffusive flux in X
# Line 208  C-    Diffusive flux in X Line 276  C-    Diffusive flux in X
276         ENDDO         ENDDO
277        ENDIF        ENDIF
278    
279    C-    Add bi-harmonic diffusive flux in X
280          IF (diffK4 .NE. 0.) THEN
281           CALL GAD_BIHARM_X(bi,bj,k,xA,df4,diffK4,df,myThid)
282          ENDIF
283    
284  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
285  C-    GM/Redi flux in X  C-    GM/Redi flux in X
286        IF (useGMRedi) THEN        IF ( trUseGMRedi ) THEN
287  C *note* should update GMREDI_XTRANSPORT to use localT and set df  *aja*  C *note* should update GMREDI_XTRANSPORT to set df  *aja*
288          CALL GMREDI_XTRANSPORT(          IF ( applyAB_onTracer ) THEN
289       I     iMin,iMax,jMin,jMax,bi,bj,K,            CALL GMREDI_XTRANSPORT(
290       I     xA,Tracer,tracerIdentity,       I         iMin,iMax,jMin,jMax,bi,bj,k,
291       U     df,       I         xA,TracerN,tracerIdentity,
292       I     myThid)       U         df,
293         I         myThid)
294            ELSE
295              CALL GMREDI_XTRANSPORT(
296         I         iMin,iMax,jMin,jMax,bi,bj,k,
297         I         xA,TracAB, tracerIdentity,
298         U         df,
299         I         myThid)
300            ENDIF
301        ENDIF        ENDIF
302  #endif  #endif
303    C     anelastic: advect.fluxes are scaled by rhoFac but hor.diff. flx are not
304        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
305         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
306          fZon(i,j) = fZon(i,j) + df(i,j)          fZon(i,j) = fZon(i,j) + df(i,j)*rhoFacC(k)
307         ENDDO         ENDDO
308        ENDDO        ENDDO
309    
310  C-    Bi-harmonic duffusive flux in X  #ifdef ALLOW_DIAGNOSTICS
311        IF (diffK4 .NE. 0.) THEN  C-    Diagnostics of Tracer flux in X dir (mainly Diffusive term),
312         CALL GAD_BIHARM_X(bi,bj,k,xA,df4,diffK4,df,myThid)  C       excluding advective terms:
313         DO j=1-Oly,sNy+Oly        IF ( useDiagnostics .AND.
314          DO i=1-Olx,sNx+Olx       &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. trUseGMRedi) ) THEN
315           fZon(i,j) = fZon(i,j) + df(i,j)            diagName = 'DFxE'//diagSufx
316          ENDDO            CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
        ENDDO  
317        ENDIF        ENDIF
318    #endif
319    
320  C--   Initialize net flux in Y direction  C--   Initialize net flux in Y direction
321        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
# Line 244  C--   Initialize net flux in Y direction Line 326  C--   Initialize net flux in Y direction
326    
327  C-    Advective flux in Y  C-    Advective flux in Y
328        IF (calcAdvection) THEN        IF (calcAdvection) THEN
329        IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN          IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
330         CALL GAD_C2_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)            CALL GAD_C2_ADV_Y(bi,bj,k,vTrans,locABT,af,myThid)
331        ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST
332         CALL GAD_FLUXLIMIT_ADV_Y(       &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
333       &       bi,bj,k,deltaTtracer,vTrans,vVel,localT,af,myThid)            CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme, .TRUE.,
334        ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN       I            dTtracerLev(k), vTrans, vFld, locABT,
335         CALL GAD_U3_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)       O            af, myThid )
336        ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
337         CALL GAD_C4_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)            CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, .TRUE., dTtracerLev(k),
338        ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN       I            vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,
339         CALL GAD_DST3_ADV_Y(       O            af, myThid )
340       &       bi,bj,k,deltaTtracer,vTrans,vVel,localT,af,myThid)          ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
341        ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN            CALL GAD_U3_ADV_Y(bi,bj,k,vTrans,locABT,af,myThid)
342         CALL GAD_DST3FL_ADV_Y(          ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
343       &       bi,bj,k,deltaTtracer,vTrans,vVel,localT,af,myThid)            CALL GAD_C4_ADV_Y(bi,bj,k,vTrans,locABT,af,myThid)
344        ELSE          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
345         STOP 'GAD_CALC_RHS: Bad advectionScheme (Y)'            CALL GAD_DST3_ADV_Y( bi,bj,k, .TRUE., dTtracerLev(k),
346        ENDIF       I            vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,
347        DO j=1-Oly,sNy+Oly       O            af, myThid )
348         DO i=1-Olx,sNx+Olx          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
349          fMer(i,j) = fMer(i,j) + af(i,j)           IF ( inAdMode ) THEN
350         ENDDO  cph This block is to trick the adjoint:
351        ENDDO  cph IF inAdExact=.FALSE., we want to use DST3
352    cph with limiters in forward, but without limiters in reverse.
353              CALL GAD_DST3_ADV_Y( bi,bj,k, .TRUE., dTtracerLev(k),
354         I           vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,
355         O           af, myThid )
356             ELSE
357              CALL GAD_DST3FL_ADV_Y( bi,bj,k, .TRUE., dTtracerLev(k),
358         I           vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,
359         O           af, myThid )
360             ENDIF
361            ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN
362              CALL GAD_OS7MP_ADV_Y( bi,bj,k, .TRUE., dTtracerLev(k),
363         I            vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,
364         O            af, myThid )
365            ELSE
366              STOP 'GAD_CALC_RHS: Bad advectionScheme (Y)'
367            ENDIF
368            DO j=1-Oly,sNy+Oly
369             DO i=1-Olx,sNx+Olx
370              fMer(i,j) = fMer(i,j) + af(i,j)
371             ENDDO
372            ENDDO
373    #ifdef ALLOW_DIAGNOSTICS
374            IF ( useDiagnostics ) THEN
375              diagName = 'ADVy'//diagSufx
376              CALL DIAGNOSTICS_FILL(af,diagName, k,1, 2,bi,bj, myThid)
377            ENDIF
378    #endif
379        ENDIF        ENDIF
380    
381  C-    Diffusive flux in Y  C-    Diffusive flux in Y
# Line 280  C-    Diffusive flux in Y Line 389  C-    Diffusive flux in Y
389         ENDDO         ENDDO
390        ENDIF        ENDIF
391    
392    C-    Add bi-harmonic flux in Y
393          IF (diffK4 .NE. 0.) THEN
394           CALL GAD_BIHARM_Y(bi,bj,k,yA,df4,diffK4,df,myThid)
395          ENDIF
396    
397  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
398  C-    GM/Redi flux in Y  C-    GM/Redi flux in Y
399        IF (useGMRedi) THEN        IF ( trUseGMRedi ) THEN
400  C *note* should update GMREDI_YTRANSPORT to use localT and set df  *aja*  C *note* should update GMREDI_YTRANSPORT to set df  *aja*
401         CALL GMREDI_YTRANSPORT(          IF ( applyAB_onTracer ) THEN
402       I     iMin,iMax,jMin,jMax,bi,bj,K,            CALL GMREDI_YTRANSPORT(
403       I     yA,Tracer,tracerIdentity,       I         iMin,iMax,jMin,jMax,bi,bj,k,
404       U     df,       I         yA,TracerN,tracerIdentity,
405       I     myThid)       U         df,
406         I         myThid)
407            ELSE
408              CALL GMREDI_YTRANSPORT(
409         I         iMin,iMax,jMin,jMax,bi,bj,k,
410         I         yA,TracAB, tracerIdentity,
411         U         df,
412         I         myThid)
413            ENDIF
414        ENDIF        ENDIF
415  #endif  #endif
416    C     anelastic: advect.fluxes are scaled by rhoFac but hor.diff. flx are not
417        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
418         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
419          fMer(i,j) = fMer(i,j) + df(i,j)          fMer(i,j) = fMer(i,j) + df(i,j)*rhoFacC(k)
420         ENDDO         ENDDO
421        ENDDO        ENDDO
422    
423  C-    Bi-harmonic flux in Y  #ifdef ALLOW_DIAGNOSTICS
424        IF (diffK4 .NE. 0.) THEN  C-    Diagnostics of Tracer flux in Y dir (mainly Diffusive terms),
425         CALL GAD_BIHARM_Y(bi,bj,k,yA,df4,diffK4,df,myThid)  C       excluding advective terms:
426         DO j=1-Oly,sNy+Oly        IF ( useDiagnostics .AND.
427          DO i=1-Olx,sNx+Olx       &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. trUseGMRedi) ) THEN
428           fMer(i,j) = fMer(i,j) + df(i,j)            diagName = 'DFyE'//diagSufx
429          ENDDO            CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
        ENDDO  
430        ENDIF        ENDIF
431    #endif
432    
433  C--   Compute vertical flux fVerT(kUp) at interface k (between k-1 & k):  C--   Compute vertical flux fVerT(kUp) at interface k (between k-1 & k):
434  C-    Advective flux in R  C-    Advective flux in R
435  #ifdef ALLOW_AIM  #ifdef ALLOW_AIM
436  C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr  C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr
437        IF (calcAdvection .AND. .NOT.implicitAdvection .AND. K.GE.2 .AND.        IF (calcAdvection .AND. .NOT.implicitAdvection .AND. k.GE.2 .AND.
438       &     (.NOT.useAIM .OR.tracerIdentity.NE.GAD_SALINITY .OR.K.LT.Nr)       &     (.NOT.useAIM .OR.tracerIdentity.NE.GAD_SALINITY .OR.k.LT.Nr)
439       &   ) THEN       &   ) THEN
440  #else  #else
441        IF (calcAdvection .AND. .NOT.implicitAdvection .AND. K.GE.2) THEN        IF (calcAdvection .AND. .NOT.implicitAdvection .AND. k.GE.2) THEN
442  #endif  #endif
443  C-    Compute vertical advective flux in the interior:  C-    Compute vertical advective flux in the interior:
444         IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN          IF (vertAdvecScheme.EQ.ENUM_CENTERED_2ND) THEN
445          CALL GAD_C2_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)            CALL GAD_C2_ADV_R(bi,bj,k,rTrans,TracAB,af,myThid)
446         ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF ( vertAdvecScheme.EQ.ENUM_UPWIND_1RST
447          CALL GAD_FLUXLIMIT_ADV_R(       &          .OR. vertAdvecScheme.EQ.ENUM_DST2 ) THEN
448       &       bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)            CALL GAD_DST2U1_ADV_R( bi,bj,k, vertAdvecScheme,
449         ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN       I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),
450          CALL GAD_U3_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)       O         af, myThid )
451         ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN          ELSEIF (vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN
452          CALL GAD_C4_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)            CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k,
453         ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN       I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),
454          CALL GAD_DST3_ADV_R(       O         af, myThid )
455       &       bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)          ELSEIF (vertAdvecScheme.EQ.ENUM_UPWIND_3RD ) THEN
456         ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN            CALL GAD_U3_ADV_R(bi,bj,k,rTrans,TracAB,af,myThid)
457          CALL GAD_DST3FL_ADV_R(          ELSEIF (vertAdvecScheme.EQ.ENUM_CENTERED_4TH) THEN
458       &       bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)            CALL GAD_C4_ADV_R(bi,bj,k,rTrans,TracAB,af,myThid)
459         ELSE          ELSEIF (vertAdvecScheme.EQ.ENUM_DST3 ) THEN
460          STOP 'GAD_CALC_RHS: Bad advectionScheme (R)'            CALL GAD_DST3_ADV_R( bi,bj,k,
461         ENDIF       I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),
462         O         af, myThid )
463            ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
464    cph This block is to trick the adjoint:
465    cph IF inAdExact=.FALSE., we want to use DST3
466    cph with limiters in forward, but without limiters in reverse.
467              IF ( inAdMode ) THEN
468               CALL GAD_DST3_ADV_R( bi,bj,k,
469         I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),
470         O         af, myThid )
471              ELSE
472               CALL GAD_DST3FL_ADV_R( bi,bj,k,
473         I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),
474         O         af, myThid )
475              ENDIF
476            ELSEIF (vertAdvecScheme.EQ.ENUM_OS7MP ) THEN
477               CALL GAD_OS7MP_ADV_R( bi,bj,k,
478         I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),
479         O         af, myThid )
480            ELSE
481              STOP 'GAD_CALC_RHS: Bad vertAdvecScheme (R)'
482            ENDIF
483  C-     add the advective flux to fVerT  C-     add the advective flux to fVerT
484         DO j=1-Oly,sNy+Oly          DO j=1-Oly,sNy+Oly
485          DO i=1-Olx,sNx+Olx           DO i=1-Olx,sNx+Olx
486           fVerT(i,j,kUp) = fVerT(i,j,kUp) + af(i,j)            fVerT(i,j,kUp) = fVerT(i,j,kUp) + af(i,j)
487             ENDDO
488          ENDDO          ENDDO
489         ENDDO  #ifdef ALLOW_DIAGNOSTICS
490            IF ( useDiagnostics ) THEN
491              diagName = 'ADVr'//diagSufx
492              CALL DIAGNOSTICS_FILL(af,diagName, k,1, 2,bi,bj, myThid)
493    C- note: needs to explicitly increment the counter since DIAGNOSTICS_FILL
494    C        does it only if k=1 (never the case here)
495              IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT(diagName,bi,bj,myThid)
496            ENDIF
497    #endif
498        ENDIF        ENDIF
499    
500  C-    Diffusive flux in R  C-    Diffusive flux in R
# Line 354  C           boundary condition. Line 507  C           boundary condition.
507          ENDDO          ENDDO
508         ENDDO         ENDDO
509        ELSE        ELSE
510         CALL GAD_DIFF_R(bi,bj,k,KappaRT,tracer,df,myThid)         IF ( applyAB_onTracer ) THEN
511             CALL GAD_DIFF_R(bi,bj,k,KappaR,TracerN,df,myThid)
512           ELSE
513             CALL GAD_DIFF_R(bi,bj,k,KappaR,TracAB, df,myThid)
514           ENDIF
515        ENDIF        ENDIF
516    
517  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
518  C-    GM/Redi flux in R  C-    GM/Redi flux in R
519        IF (useGMRedi) THEN        IF ( trUseGMRedi ) THEN
520  C *note* should update GMREDI_RTRANSPORT to set df  *aja*  C *note* should update GMREDI_RTRANSPORT to set df  *aja*
521         CALL GMREDI_RTRANSPORT(          IF ( applyAB_onTracer ) THEN
522       I     iMin,iMax,jMin,jMax,bi,bj,K,            CALL GMREDI_RTRANSPORT(
523       I     Tracer,tracerIdentity,       I         iMin,iMax,jMin,jMax,bi,bj,k,
524       U     df,       I         TracerN,tracerIdentity,
525       I     myThid)       U         df,
526         I         myThid)
527            ELSE
528              CALL GMREDI_RTRANSPORT(
529         I         iMin,iMax,jMin,jMax,bi,bj,k,
530         I         TracAB, tracerIdentity,
531         U         df,
532         I         myThid)
533            ENDIF
534        ENDIF        ENDIF
535  #endif  #endif
536    
# Line 375  C *note* should update GMREDI_RTRANSPORT Line 540  C *note* should update GMREDI_RTRANSPORT
540         ENDDO         ENDDO
541        ENDDO        ENDDO
542    
543    #ifdef ALLOW_DIAGNOSTICS
544    C-    Diagnostics of Tracer flux in R dir (mainly Diffusive terms),
545    C       Explicit terms only & excluding advective terms:
546          IF ( useDiagnostics .AND.
547         &    (.NOT.implicitDiffusion .OR. trUseGMRedi) ) THEN
548              diagName = 'DFrE'//diagSufx
549              CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
550          ENDIF
551    #endif
552    
553  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
554  C-    Add non local KPP transport term (ghat) to diffusive T flux.  C-    Set non local KPP transport term (ghat):
555        IF (useKPP) THEN        IF ( trUseKPP .AND. k.GE.2 ) THEN
556         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
557          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
558           df(i,j) = 0. _d 0           df(i,j) = 0. _d 0
559          ENDDO          ENDDO
560         ENDDO         ENDDO
561         IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN         IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
 C *note* should update KPP_TRANSPORT_T to set df  *aja*  
562          CALL KPP_TRANSPORT_T(          CALL KPP_TRANSPORT_T(
563       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,       I           iMin,iMax,jMin,jMax,bi,bj,k,km1,
564       I     KappaRT,       O           df,
565       U     df )       I           myTime, myIter, myThid )
566         ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN         ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
567          CALL KPP_TRANSPORT_S(          CALL KPP_TRANSPORT_S(
568       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,       I           iMin,iMax,jMin,jMax,bi,bj,k,km1,
569       I     KappaRT,       O           df,
570       U     df )       I           myTime, myIter, myThid )
571  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
572         ELSEIF (tracerIdentity .GE. GAD_TR1) THEN         ELSEIF (tracerIdentity .GE. GAD_TR1) THEN
573          CALL KPP_TRANSPORT_PTR(          CALL KPP_TRANSPORT_PTR(
574       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,       I           iMin,iMax,jMin,jMax,bi,bj,k,km1,
575       I     tracerIdentity-GAD_TR1+1,KappaRT,       I           tracerIdentity-GAD_TR1+1,
576       U     df )       O           df,
577         I           myTime, myIter, myThid )
578  #endif  #endif
579         ELSE         ELSE
580          PRINT*,'invalid tracer indentity: ', tracerIdentity          PRINT*,'invalid tracer indentity: ', tracerIdentity
# Line 407  C *note* should update KPP_TRANSPORT_T t Line 582  C *note* should update KPP_TRANSPORT_T t
582         ENDIF         ENDIF
583         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
584          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
585           fVerT(i,j,kUp) = fVerT(i,j,kUp) + df(i,j)*maskUp(i,j)           fVerT(i,j,kUp) = fVerT(i,j,kUp)
586         &                  + df(i,j)*maskUp(i,j)*rhoFacF(k)
587            ENDDO
588           ENDDO
589          ENDIF
590    #endif
591    
592    #ifdef GAD_SMOLARKIEWICZ_HACK
593    coj   Hack to make redi (and everything else in this s/r) positive
594    coj   (see Smolarkiewicz MWR 1989 and Bott MWR 1989).
595    coj   Only works if 'down' is k+1 and k loop in thermodynamics is k=Nr,1,-1
596    coj
597    coj   Apply to all tracers except temperature
598          IF (tracerIdentity.NE.GAD_TEMPERATURE .AND.
599         &    tracerIdentity.NE.GAD_SALINITY) THEN
600           DO j=1-Oly,sNy+Oly-1
601            DO i=1-Olx,sNx+Olx-1
602    coj   Add outgoing fluxes
603             outFlux=dTtracerLev(k)*
604         &    _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
605         &   *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)*recip_rhoFacC(k)
606         &    *( MAX(0. _d 0,fZon(i+1,j)) + MAX(0. _d 0,-fZon(i,j))
607         &      +MAX(0. _d 0,fMer(i,j+1)) + MAX(0. _d 0,-fMer(i,j))
608         &      +MAX(0. _d 0,fVerT(i,j,kDown)*rkSign)
609         &      +MAX(0. _d 0,-fVerT(i,j,kUp)*rkSign)
610         &     )
611             IF ( applyAB_onTracer ) THEN
612               trac=TracerN(i,j,k,bi,bj)
613             ELSE
614               trac=TracAB(i,j,k,bi,bj)
615             ENDIF
616    coj   If they would reduce tracer by a fraction of more than
617    coj   SmolarkiewiczMaxFrac, scale them down
618             IF (outFlux.GT.0. _d 0 .AND.
619         &       outFlux.GT.SmolarkiewiczMaxFrac*trac) THEN
620    coj   If tracer is already negative, scale flux to zero
621               fac = MAX(0. _d 0,SmolarkiewiczMaxFrac*trac/outFlux)
622    
623               IF (fZon(i+1,j).GT.0. _d 0) fZon(i+1,j)=fac*fZon(i+1,j)
624               IF (-fZon(i,j) .GT.0. _d 0) fZon(i,j)  =fac*fZon(i,j)      
625               IF (fMer(i,j+1).GT.0. _d 0) fMer(i,j+1)=fac*fMer(i,j+1)
626               IF (-fMer(i,j) .GT.0. _d 0) fMer(i,j)  =fac*fMer(i,j)
627               IF (-fVerT(i,j,kUp)*rkSign .GT.0. _d 0)
628         &       fVerT(i,j,kUp)=fac*fVerT(i,j,kUp)
629    
630               IF (k.LT.Nr .AND. fVerT(i,j,kDown)*rkSign.GT.0. _d 0) THEN
631    coj   Down flux is special: it has already been applied in lower layer,
632    coj   so we have to readjust this.
633    coj   Note: for k+1, gTracer is now the updated tracer, not the tendency!
634    coj   thus it has an extra factor dTtracerLev(k+1)
635                 gTrFac=dTtracerLev(k+1)
636    coj   Other factors that have been applied to gTracer since the last call:
637    #ifdef NONLIN_FRSURF
638                 IF (nonlinFreeSurf.GT.0) THEN
639                  IF (select_rStar.GT.0) THEN
640    #ifndef DISABLE_RSTAR_CODE
641                    gTrFac = gTrFac/rStarExpC(i,j,bi,bj)
642    #endif /* DISABLE_RSTAR_CODE */
643                  ENDIF
644                 ENDIF
645    #endif /* NONLIN_FRSURF */
646    coj   Now: undo down flux, ...
647                 gTracer(i,j,k+1,bi,bj)=gTracer(i,j,k+1,bi,bj)
648         &        +gTrFac
649         &         *_recip_hFacC(i,j,k+1,bi,bj)*recip_drF(k+1)
650         &         *recip_rA(i,j,bi,bj)*recip_deepFac2C(k+1)
651         &         *recip_rhoFacC(k+1)
652         &         *( -fVerT(i,j,kDown)*rkSign )
653    coj   ... scale ...
654                 fVerT(i,j,kDown)=fac*fVerT(i,j,kDown)
655    coj   ... and reapply
656                 gTracer(i,j,k+1,bi,bj)=gTracer(i,j,k+1,bi,bj)
657         &        +gTrFac
658         &         *_recip_hFacC(i,j,k+1,bi,bj)*recip_drF(k+1)
659         &         *recip_rA(i,j,bi,bj)*recip_deepFac2C(k+1)
660         &         *recip_rhoFacC(k+1)
661         &         *( fVerT(i,j,kDown)*rkSign )
662               ENDIF
663    
664             ENDIF
665          ENDDO          ENDDO
666         ENDDO         ENDDO
667        ENDIF        ENDIF
668  #endif  #endif
669    
670  C--   Divergence of fluxes  C--   Divergence of fluxes
671    C     Anelastic: scale vertical fluxes by rhoFac and leave Horizontal fluxes unchanged
672        DO j=1-Oly,sNy+Oly-1        DO j=1-Oly,sNy+Oly-1
673         DO i=1-Olx,sNx+Olx-1         DO i=1-Olx,sNx+Olx-1
674          gTracer(i,j,k,bi,bj)=gTracer(i,j,k,bi,bj)          gTracer(i,j,k,bi,bj)=gTracer(i,j,k,bi,bj)
675       &   -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)*recip_rA(i,j,bi,bj)       &   -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
676         &   *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)*recip_rhoFacC(k)
677       &   *( (fZon(i+1,j)-fZon(i,j))       &   *( (fZon(i+1,j)-fZon(i,j))
678       &     +(fMer(i,j+1)-fMer(i,j))       &     +(fMer(i,j+1)-fMer(i,j))
679       &     +(fVerT(i,j,kUp)-fVerT(i,j,kDown))*rkFac       &     +(fVerT(i,j,kDown)-fVerT(i,j,kUp))*rkSign
680       &     -localT(i,j)*( (uTrans(i+1,j)-uTrans(i,j))       &     -localT(i,j)*( (uTrans(i+1,j)-uTrans(i,j))
681       &                   +(vTrans(i,j+1)-vTrans(i,j))       &                   +(vTrans(i,j+1)-vTrans(i,j))
682       &                   +(rTrans(i,j)-rTransKp1(i,j))*rAdvFac       &                   +(rTransKp1(i,j)-rTrans(i,j))*rAdvFac
683       &                  )*advFac       &                  )*advFac
684       &    )       &    )
685         ENDDO         ENDDO
686        ENDDO        ENDDO
687    
688    #ifdef ALLOW_DEBUG
689          IF ( debugLevel .GE. debLevB
690         &   .AND. tracerIdentity.EQ.GAD_TEMPERATURE
691         &   .AND. k.EQ.2 .AND. myIter.EQ.1+nIter0
692         &   .AND. nPx.EQ.1 .AND. nPy.EQ.1
693         &   .AND. useCubedSphereExchange ) THEN
694            CALL DEBUG_CS_CORNER_UV( ' fZon,fMer from GAD_CALC_RHS',
695         &             fZon,fMer, k, standardMessageUnit,bi,bj,myThid )
696          ENDIF
697    #endif /* ALLOW_DEBUG */
698    
699        RETURN        RETURN
700        END        END
Legend:
Removed from v.1.25  
changed lines
  Added in v.1.52
  ViewVC Help
Powered by ViewVC 1.1.22