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