/[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.17 by jmc, Sun Jan 26 21:08:36 2003 UTC revision 1.61 by jmc, Fri Jun 15 13:28:41 2012 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,maskUp,       I           xA, yA, maskUp, uFld, vFld, wFld,
13       I           diffKh, diffK4, KappaRT, Tracer,       I           uTrans, vTrans, rTrans, rTransKp1,
14       I           tracerIdentity, advectionScheme, calcAdvection,       I           diffKh, diffK4, KappaR, diffKr4, TracerN, TracAB,
15         I           deltaTLev, trIdentity,
16         I           advectionScheme, vertAdvecScheme,
17         I           calcAdvection, implicitAdvection, applyAB_onTracer,
18         I           trUseDiffKr4, trUseGMRedi, trUseKPP,
19       U           fVerT, gTracer,       U           fVerT, gTracer,
20       I           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 27  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 40  C !USES: =============================== Line 44  C !USES: ===============================
44  #include "EEPARAMS.h"  #include "EEPARAMS.h"
45  #include "PARAMS.h"  #include "PARAMS.h"
46  #include "GRID.h"  #include "GRID.h"
 #include "DYNVARS.h"  
47  #include "SURFACE.h"  #include "SURFACE.h"
48  #include "GAD.h"  #include "GAD.h"
49    
# Line 50  C !USES: =============================== Line 53  C !USES: ===============================
53  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
54    
55  C !INPUT PARAMETERS: ===================================================  C !INPUT PARAMETERS: ===================================================
56  C  bi,bj                :: tile indices  C bi,bj            :: tile indices
57  C  iMin,iMax,jMin,jMax  :: loop range for called routines  C iMin,iMax        :: loop range for called routines
58  C  kup                  :: index into 2 1/2D array, toggles between 1 and 2  C jMin,jMax        :: loop range for called routines
59  C  kdown                :: index into 2 1/2D array, toggles between 2 and 1  C k                :: vertical index
60  C  kp1                  :: =k+1 for k<Nr, =Nr for k=Nr  C kM1              :: =k-1 for k>1, =1 for k=1
61  C  xA,yA                :: areas of X and Y face of tracer cells  C kUp              :: index into 2 1/2D array, toggles between 1|2
62  C  uTrans,vTrans,rTrans :: 2-D arrays of volume transports at U,V and W points  C kDown            :: index into 2 1/2D array, toggles between 2|1
63  C  maskUp               :: 2-D array for mask at W points  C xA,yA            :: areas of X and Y face of tracer cells
64  C  diffKh               :: horizontal diffusion coefficient  C maskUp           :: 2-D array for mask at W points
65  C  diffK4               :: bi-harmonic diffusion coefficient  C uFld,vFld,wFld   :: Local copy of velocity field (3 components)
66  C  KappaRT              :: 3-D array for vertical diffusion coefficient  C uTrans,vTrans    :: 2-D arrays of volume transports at U,V points
67  C  Tracer               :: tracer field  C rTrans           :: 2-D arrays of volume transports at W points
68  C  tracerIdentity       :: identifier for the tracer (required only for KPP)  C rTransKp1        :: 2-D array of volume trans at W pts, interf k+1
69  C  advectionScheme      :: advection scheme to use  C diffKh           :: horizontal diffusion coefficient
70  C  calcAdvection        :: =False if Advec terms computed with multiDim scheme  C diffK4           :: horizontal bi-harmonic diffusion coefficient
71  C  myThid               :: thread number  C KappaR           :: 2-D array for vertical diffusion coefficient, interf k
72    C diffKr4          :: 1-D array for vertical bi-harmonic diffusion coefficient
73    C TracerN          :: tracer field @ time-step n (Note: only used
74    C                     if applying AB on tracer field rather than on tendency gTr)
75    C TracAB           :: current tracer field (@ time-step n if applying AB on gTr
76    C                     or extrapolated fwd in time to n+1/2 if applying AB on Tr)
77    C trIdentity       :: tracer identifier (required for KPP,GM)
78    C advectionScheme  :: advection scheme to use (Horizontal plane)
79    C vertAdvecScheme  :: advection scheme to use (Vertical direction)
80    C calcAdvection    :: =False if Advec computed with multiDim scheme
81    C implicitAdvection:: =True if vertical Advec computed implicitly
82    C applyAB_onTracer :: apply Adams-Bashforth on Tracer (rather than on gTr)
83    C trUseDiffKr4     :: true if this tracer uses vertical bi-harmonic diffusion
84    C trUseGMRedi      :: true if this tracer uses GM-Redi
85    C trUseKPP         :: true if this tracer uses KPP
86    C myTime           :: current time
87    C myIter           :: iteration number
88    C myThid           :: thread number
89        INTEGER bi,bj,iMin,iMax,jMin,jMax        INTEGER bi,bj,iMin,iMax,jMin,jMax
90        INTEGER k,kUp,kDown,kM1        INTEGER k,kUp,kDown,kM1
91        _RS xA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS xA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
92        _RS yA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS yA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
93          _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
94          _RL uFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
95          _RL vFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
96          _RL wFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
97        _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
98        _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
99        _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
100        _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
101        _RL diffKh, diffK4        _RL diffKh, diffK4
102        _RL KappaRT(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL KappaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
103        _RL Tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL diffKr4(Nr)
104        INTEGER tracerIdentity        _RL TracerN(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
105        INTEGER advectionScheme        _RL TracAB (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
106          _RL deltaTLev(Nr)
107          INTEGER trIdentity
108          INTEGER advectionScheme, vertAdvecScheme
109        LOGICAL calcAdvection        LOGICAL calcAdvection
110        INTEGER myThid        LOGICAL implicitAdvection, applyAB_onTracer
111          LOGICAL trUseDiffKr4, trUseGMRedi, trUseKPP
112          _RL     myTime
113          INTEGER myIter, myThid
114    
115  C !OUTPUT PARAMETERS: ==================================================  C !OUTPUT PARAMETERS: ==================================================
116  C  gTracer              :: tendancy array  C gTracer          :: tendency array
117  C  fVerT                :: 2 1/2D arrays for vertical advective flux  C fVerT            :: 2 1/2D arrays for vertical advective flux
118        _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)
119        _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
120    
121    C !FUNCTIONS:       ====================================================
122    #ifdef ALLOW_DIAGNOSTICS
123          CHARACTER*4 GAD_DIAG_SUFX
124          EXTERNAL    GAD_DIAG_SUFX
125    #endif /* ALLOW_DIAGNOSTICS */
126    
127  C !LOCAL VARIABLES: ====================================================  C !LOCAL VARIABLES: ====================================================
128  C  i,j                  :: loop indices  C i,j              :: loop indices
129  C  df4                  :: used for storing del^2 T for bi-harmonic term  C df4              :: used for storing del^2 T for bi-harmonic term
130  C  fZon                 :: zonal flux  C fZon             :: zonal flux
131  C  fmer                 :: meridional flux  C fMer             :: meridional flux
132  C  af                   :: advective flux  C af               :: advective flux
133  C  df                   :: diffusive flux  C df               :: diffusive flux
134  C  localT               :: local copy of tracer field  C localT           :: local copy of tracer field
135    C locABT           :: local copy of (AB-extrapolated) tracer field
136        INTEGER i,j        INTEGER i,j
137          _RS maskLocW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
138          _RS maskLocS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
139        _RL df4   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL df4   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
140        _RL fZon  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fZon  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
141        _RL fMer  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fMer  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
142        _RL af    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL af    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
143        _RL df    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL df    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
144        _RL localT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL localT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
145          _RL locABT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
146          _RL advFac, rAdvFac
147    #ifdef GAD_SMOLARKIEWICZ_HACK
148          _RL outFlux, trac, fac, gTrFac
149    #endif
150    #ifdef ALLOW_DIAGNOSTICS
151          CHARACTER*8 diagName
152          CHARACTER*4 diagSufx
153    #endif
154  CEOP  CEOP
155    
156  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 111  C--   the kDown is still required Line 159  C--   the kDown is still required
159        fVerT(1,1,kDown) = fVerT(1,1,kDown)        fVerT(1,1,kDown) = fVerT(1,1,kDown)
160  #endif  #endif
161    
162    #ifdef ALLOW_DIAGNOSTICS
163    C--   Set diagnostic suffix for the current tracer
164          IF ( useDiagnostics ) THEN
165            diagSufx = GAD_DIAG_SUFX( trIdentity, myThid )
166          ENDIF
167    #endif
168    
169          advFac  = 0. _d 0
170          IF (calcAdvection) advFac = 1. _d 0
171          rAdvFac = rkSign*advFac
172          IF (implicitAdvection) rAdvFac = rkSign
173    
174        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
175         DO i=1-OLx,sNx+OLx         DO i=1-OLx,sNx+OLx
176          fZon(i,j)      = 0. _d 0          fZon(i,j)      = 0. _d 0
# Line 118  C--   the kDown is still required Line 178  C--   the kDown is still required
178          fVerT(i,j,kUp) = 0. _d 0          fVerT(i,j,kUp) = 0. _d 0
179          df(i,j)        = 0. _d 0          df(i,j)        = 0. _d 0
180          df4(i,j)       = 0. _d 0          df4(i,j)       = 0. _d 0
         localT(i,j)    = 0. _d 0  
181         ENDDO         ENDDO
182        ENDDO        ENDDO
183    
184  C--   Make local copy of tracer array  C--   Make local copy of tracer array
185        DO j=1-OLy,sNy+OLy        IF ( applyAB_onTracer ) THEN
186         DO i=1-OLx,sNx+OLx          DO j=1-OLy,sNy+OLy
187          localT(i,j)=tracer(i,j,k,bi,bj)           DO i=1-OLx,sNx+OLx
188         ENDDO            localT(i,j)=TracerN(i,j,k,bi,bj)
189        ENDDO            locABT(i,j)= TracAB(i,j,k,bi,bj)
190             ENDDO
191            ENDDO
192          ELSE
193            DO j=1-OLy,sNy+OLy
194             DO i=1-OLx,sNx+OLx
195              localT(i,j)= TracAB(i,j,k,bi,bj)
196              locABT(i,j)= TracAB(i,j,k,bi,bj)
197             ENDDO
198            ENDDO
199          ENDIF
200    
201  C--   Unless we have already calculated the advection terms we initialize  C--   Unless we have already calculated the advection terms we initialize
202  C     the tendency to zero.  C     the tendency to zero.
203        IF (calcAdvection) THEN  C     <== now done earlier at the beginning of thermodynamics.
204         DO j=1-Oly,sNy+Oly  c     IF (calcAdvection) THEN
205          DO i=1-Olx,sNx+Olx  c      DO j=1-OLy,sNy+OLy
206           gTracer(i,j,k,bi,bj)=0. _d 0  c       DO i=1-OLx,sNx+OLx
207          ENDDO  c        gTracer(i,j,k,bi,bj)=0. _d 0
208         ENDDO  c       ENDDO
209        ENDIF  c      ENDDO
210    c     ENDIF
211    
212  C--   Pre-calculate del^2 T if bi-harmonic coefficient is non-zero  C--   Pre-calculate del^2 T if bi-harmonic coefficient is non-zero
213        IF (diffK4 .NE. 0.) THEN        IF (diffK4 .NE. 0.) THEN
# Line 147  C--   Pre-calculate del^2 T if bi-harmon Line 217  C--   Pre-calculate del^2 T if bi-harmon
217        ENDIF        ENDIF
218    
219  C--   Initialize net flux in X direction  C--   Initialize net flux in X direction
220        DO j=1-Oly,sNy+Oly        DO j=1-OLy,sNy+OLy
221         DO i=1-Olx,sNx+Olx         DO i=1-OLx,sNx+OLx
222          fZon(i,j) = 0. _d 0          fZon(i,j) = 0. _d 0
223         ENDDO         ENDDO
224        ENDDO        ENDDO
225    
226  C-    Advective flux in X  C-    Advective flux in X
227        IF (calcAdvection) THEN        IF (calcAdvection) THEN
228        IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN          IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
229         CALL GAD_C2_ADV_X(bi,bj,k,uTrans,localT,af,myThid)             CALL GAD_C2_ADV_X( bi,bj,k, uTrans, locABT, af, myThid )
230        ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST
231         CALL GAD_FLUXLIMIT_ADV_X(       &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
232       &      bi,bj,k,deltaTtracer,uTrans,uVel,localT,af,myThid)             CALL GAD_DST2U1_ADV_X( bi,bj,k, advectionScheme, .TRUE.,
233        ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN       I              deltaTLev(k), uTrans, uFld, locABT,
234         CALL GAD_U3_ADV_X(bi,bj,k,uTrans,localT,af,myThid)       O              af, myThid )
235        ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN          ELSE
236         CALL GAD_C4_ADV_X(bi,bj,k,uTrans,localT,af,myThid)           DO j=1-OLy,sNy+OLy
237        ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN            DO i=1-OLx,sNx+OLx
238         CALL GAD_DST3_ADV_X(  #ifdef ALLOW_OBCS
239       &       bi,bj,k,deltaTtracer,uTrans,uVel,localT,af,myThid)             maskLocW(i,j) = _maskW(i,j,k,bi,bj)*maskInW(i,j,bi,bj)
240        ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN  #else /* ALLOW_OBCS */
241         CALL GAD_DST3FL_ADV_X(             maskLocW(i,j) = _maskW(i,j,k,bi,bj)
242       &       bi,bj,k,deltaTtracer,uTrans,uVel,localT,af,myThid)  #endif /* ALLOW_OBCS */
243        ELSE            ENDDO
244         STOP 'GAD_CALC_RHS: Bad advectionScheme (X)'           ENDDO
245        ENDIF           IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
246        DO j=1-Oly,sNy+Oly             CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
247         DO i=1-Olx,sNx+Olx       I              uTrans, uFld, maskLocW, locABT,
248          fZon(i,j) = fZon(i,j) + af(i,j)       O              af, myThid )
249         ENDDO           ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
250        ENDDO             CALL GAD_U3_ADV_X( bi,bj,k, uTrans, maskLocW, locABT,
251         O              af, myThid )
252             ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
253               CALL GAD_C4_ADV_X( bi,bj,k, uTrans, maskLocW, locABT,
254         O              af, myThid )
255             ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
256               CALL GAD_DST3_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
257         I              uTrans, uFld, maskLocW, locABT,
258         O              af, myThid )
259             ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
260              IF ( inAdMode ) THEN
261    cph This block is to trick the adjoint:
262    cph IF inAdExact=.FALSE., we want to use DST3
263    cph with limiters in forward, but without limiters in reverse.
264               CALL GAD_DST3_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
265         I              uTrans, uFld, maskLocW, locABT,
266         O              af, myThid )
267              ELSE
268               CALL GAD_DST3FL_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
269         I              uTrans, uFld, maskLocW, locABT,
270         O              af, myThid )
271              ENDIF
272    #ifndef ALLOW_AUTODIFF_TAMC
273             ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN
274               CALL GAD_OS7MP_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
275         I              uTrans, uFld, maskLocW, locABT,
276         O              af, myThid )
277    #endif
278             ELSE
279              STOP 'GAD_CALC_RHS: Bad advectionScheme (X)'
280             ENDIF
281            ENDIF
282    #ifdef ALLOW_OBCS
283            IF ( useOBCS ) THEN
284    C-      replace advective flux with 1st order upwind scheme estimate
285              CALL OBCS_U1_ADV_TRACER( .TRUE., trIdentity, bi, bj, k,
286         I                             maskW(1-OLx,1-OLy,k,bi,bj),
287         I                             uTrans, locABT,
288         U                             af, myThid )
289            ENDIF
290    #endif /* ALLOW_OBCS */
291            DO j=1-OLy,sNy+OLy
292             DO i=1-OLx,sNx+OLx
293              fZon(i,j) = fZon(i,j) + af(i,j)
294             ENDDO
295            ENDDO
296    #ifdef ALLOW_DIAGNOSTICS
297            IF ( useDiagnostics ) THEN
298              diagName = 'ADVx'//diagSufx
299              CALL DIAGNOSTICS_FILL( af, diagName, k,1, 2,bi,bj, myThid )
300            ENDIF
301    #endif
302        ENDIF        ENDIF
303    
304  C-    Diffusive flux in X  C-    Diffusive flux in X
305        IF (diffKh.NE.0.) THEN        IF (diffKh.NE.0.) THEN
306         CALL GAD_DIFF_X(bi,bj,k,xA,diffKh,localT,df,myThid)         CALL GAD_DIFF_X(bi,bj,k,xA,diffKh,localT,df,myThid)
307        ELSE        ELSE
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           df(i,j) = 0. _d 0           df(i,j) = 0. _d 0
311          ENDDO          ENDDO
312         ENDDO         ENDDO
313        ENDIF        ENDIF
314    
315    C-    Add bi-harmonic diffusive flux in X
316          IF (diffK4 .NE. 0.) THEN
317           CALL GAD_BIHARM_X(bi,bj,k,xA,df4,diffK4,df,myThid)
318          ENDIF
319    
320  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
321  C-    GM/Redi flux in X  C-    GM/Redi flux in X
322        IF (useGMRedi) THEN        IF ( trUseGMRedi ) THEN
323  C *note* should update GMREDI_XTRANSPORT to use localT and set df  *aja*  C *note* should update GMREDI_XTRANSPORT to set df  *aja*
324          CALL GMREDI_XTRANSPORT(          IF ( applyAB_onTracer ) THEN
325       I     iMin,iMax,jMin,jMax,bi,bj,K,            CALL GMREDI_XTRANSPORT(
326       I     xA,Tracer,tracerIdentity,       I         iMin,iMax,jMin,jMax,bi,bj,k,
327       U     df,       I         xA,TracerN,trIdentity,
328       I     myThid)       U         df,
329         I         myThid)
330            ELSE
331              CALL GMREDI_XTRANSPORT(
332         I         iMin,iMax,jMin,jMax,bi,bj,k,
333         I         xA,TracAB, trIdentity,
334         U         df,
335         I         myThid)
336            ENDIF
337        ENDIF        ENDIF
338  #endif  #endif
339        DO j=1-Oly,sNy+Oly  C     anelastic: advect.fluxes are scaled by rhoFac but hor.diff. flx are not
340         DO i=1-Olx,sNx+Olx        DO j=1-OLy,sNy+OLy
341          fZon(i,j) = fZon(i,j) + df(i,j)         DO i=1-OLx,sNx+OLx
342            fZon(i,j) = fZon(i,j) + df(i,j)*rhoFacC(k)
343         ENDDO         ENDDO
344        ENDDO        ENDDO
345    
346  C-    Bi-harmonic duffusive flux in X  #ifdef ALLOW_DIAGNOSTICS
347        IF (diffK4 .NE. 0.) THEN  C-    Diagnostics of Tracer flux in X dir (mainly Diffusive term),
348         CALL GAD_BIHARM_X(bi,bj,k,xA,df4,diffK4,df,myThid)  C       excluding advective terms:
349         DO j=1-Oly,sNy+Oly        IF ( useDiagnostics .AND.
350          DO i=1-Olx,sNx+Olx       &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. trUseGMRedi) ) THEN
351           fZon(i,j) = fZon(i,j) + df(i,j)            diagName = 'DFxE'//diagSufx
352          ENDDO            CALL DIAGNOSTICS_FILL( df, diagName, k,1, 2,bi,bj, myThid )
        ENDDO  
353        ENDIF        ENDIF
354    #endif
355    
356  C--   Initialize net flux in Y direction  C--   Initialize net flux in Y direction
357        DO j=1-Oly,sNy+Oly        DO j=1-OLy,sNy+OLy
358         DO i=1-Olx,sNx+Olx         DO i=1-OLx,sNx+OLx
359          fMer(i,j) = 0. _d 0          fMer(i,j) = 0. _d 0
360         ENDDO         ENDDO
361        ENDDO        ENDDO
362    
363  C-    Advective flux in Y  C-    Advective flux in Y
364        IF (calcAdvection) THEN        IF (calcAdvection) THEN
365        IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN          IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
366         CALL GAD_C2_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)             CALL GAD_C2_ADV_Y( bi,bj,k, vTrans, locABT, af, myThid )
367        ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST
368         CALL GAD_FLUXLIMIT_ADV_Y(       &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
369       &       bi,bj,k,deltaTtracer,vTrans,vVel,localT,af,myThid)             CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme, .TRUE.,
370        ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN       I              deltaTLev(k), vTrans, vFld, locABT,
371         CALL GAD_U3_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)       O              af, myThid )
372        ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN          ELSE
373         CALL GAD_C4_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)           DO j=1-OLy,sNy+OLy
374        ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN            DO i=1-OLx,sNx+OLx
375         CALL GAD_DST3_ADV_Y(  #ifdef ALLOW_OBCS
376       &       bi,bj,k,deltaTtracer,vTrans,vVel,localT,af,myThid)             maskLocS(i,j) = _maskS(i,j,k,bi,bj)*maskInS(i,j,bi,bj)
377        ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN  #else /* ALLOW_OBCS */
378         CALL GAD_DST3FL_ADV_Y(             maskLocS(i,j) = _maskS(i,j,k,bi,bj)
379       &       bi,bj,k,deltaTtracer,vTrans,vVel,localT,af,myThid)  #endif /* ALLOW_OBCS */
380        ELSE            ENDDO
381         STOP 'GAD_CALC_RHS: Bad advectionScheme (Y)'           ENDDO
382        ENDIF           IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
383        DO j=1-Oly,sNy+Oly             CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
384         DO i=1-Olx,sNx+Olx       I              vTrans, vFld, maskLocS, locABT,
385          fMer(i,j) = fMer(i,j) + af(i,j)       O              af, myThid )
386         ENDDO           ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
387        ENDDO             CALL GAD_U3_ADV_Y( bi,bj,k, vTrans, maskLocS, locABT,
388         O              af, myThid )
389             ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
390               CALL GAD_C4_ADV_Y( bi,bj,k, vTrans, maskLocS, locABT,
391         O              af, myThid )
392             ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
393               CALL GAD_DST3_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
394         I              vTrans, vFld, maskLocS, locABT,
395         O              af, myThid )
396             ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
397              IF ( inAdMode ) THEN
398    cph This block is to trick the adjoint:
399    cph IF inAdExact=.FALSE., we want to use DST3
400    cph with limiters in forward, but without limiters in reverse.
401               CALL GAD_DST3_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
402         I              vTrans, vFld, maskLocS, locABT,
403         O              af, myThid )
404              ELSE
405               CALL GAD_DST3FL_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
406         I              vTrans, vFld, maskLocS, locABT,
407         O              af, myThid )
408              ENDIF
409    #ifndef ALLOW_AUTODIFF_TAMC
410             ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN
411               CALL GAD_OS7MP_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
412         I              vTrans, vFld, maskLocS, locABT,
413         O              af, myThid )
414    #endif
415             ELSE
416               STOP 'GAD_CALC_RHS: Bad advectionScheme (Y)'
417             ENDIF
418            ENDIF
419    #ifdef ALLOW_OBCS
420            IF ( useOBCS ) THEN
421    C-      replace advective flux with 1st order upwind scheme estimate
422              CALL OBCS_U1_ADV_TRACER( .FALSE., trIdentity, bi, bj, k,
423         I                             maskS(1-OLx,1-OLy,k,bi,bj),
424         I                             vTrans, locABT,
425         U                             af, myThid )
426            ENDIF
427    #endif /* ALLOW_OBCS */
428            DO j=1-OLy,sNy+OLy
429             DO i=1-OLx,sNx+OLx
430              fMer(i,j) = fMer(i,j) + af(i,j)
431             ENDDO
432            ENDDO
433    #ifdef ALLOW_DIAGNOSTICS
434            IF ( useDiagnostics ) THEN
435              diagName = 'ADVy'//diagSufx
436              CALL DIAGNOSTICS_FILL( af, diagName, k,1, 2,bi,bj, myThid )
437            ENDIF
438    #endif
439        ENDIF        ENDIF
440    
441  C-    Diffusive flux in Y  C-    Diffusive flux in Y
442        IF (diffKh.NE.0.) THEN        IF (diffKh.NE.0.) THEN
443         CALL GAD_DIFF_Y(bi,bj,k,yA,diffKh,localT,df,myThid)         CALL GAD_DIFF_Y(bi,bj,k,yA,diffKh,localT,df,myThid)
444        ELSE        ELSE
445         DO j=1-Oly,sNy+Oly         DO j=1-OLy,sNy+OLy
446          DO i=1-Olx,sNx+Olx          DO i=1-OLx,sNx+OLx
447           df(i,j) = 0. _d 0           df(i,j) = 0. _d 0
448          ENDDO          ENDDO
449         ENDDO         ENDDO
450        ENDIF        ENDIF
451    
452    C-    Add bi-harmonic flux in Y
453          IF (diffK4 .NE. 0.) THEN
454           CALL GAD_BIHARM_Y(bi,bj,k,yA,df4,diffK4,df,myThid)
455          ENDIF
456    
457  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
458  C-    GM/Redi flux in Y  C-    GM/Redi flux in Y
459        IF (useGMRedi) THEN        IF ( trUseGMRedi ) THEN
460  C *note* should update GMREDI_YTRANSPORT to use localT and set df  *aja*  C *note* should update GMREDI_YTRANSPORT to set df  *aja*
461         CALL GMREDI_YTRANSPORT(          IF ( applyAB_onTracer ) THEN
462       I     iMin,iMax,jMin,jMax,bi,bj,K,            CALL GMREDI_YTRANSPORT(
463       I     yA,Tracer,tracerIdentity,       I         iMin,iMax,jMin,jMax,bi,bj,k,
464       U     df,       I         yA,TracerN,trIdentity,
465       I     myThid)       U         df,
466         I         myThid)
467            ELSE
468              CALL GMREDI_YTRANSPORT(
469         I         iMin,iMax,jMin,jMax,bi,bj,k,
470         I         yA,TracAB, trIdentity,
471         U         df,
472         I         myThid)
473            ENDIF
474        ENDIF        ENDIF
475  #endif  #endif
476        DO j=1-Oly,sNy+Oly  C     anelastic: advect.fluxes are scaled by rhoFac but hor.diff. flx are not
477         DO i=1-Olx,sNx+Olx        DO j=1-OLy,sNy+OLy
478          fMer(i,j) = fMer(i,j) + df(i,j)         DO i=1-OLx,sNx+OLx
479            fMer(i,j) = fMer(i,j) + df(i,j)*rhoFacC(k)
480         ENDDO         ENDDO
481        ENDDO        ENDDO
482    
483  C-    Bi-harmonic flux in Y  #ifdef ALLOW_DIAGNOSTICS
484        IF (diffK4 .NE. 0.) THEN  C-    Diagnostics of Tracer flux in Y dir (mainly Diffusive terms),
485         CALL GAD_BIHARM_Y(bi,bj,k,yA,df4,diffK4,df,myThid)  C       excluding advective terms:
486         DO j=1-Oly,sNy+Oly        IF ( useDiagnostics .AND.
487          DO i=1-Olx,sNx+Olx       &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. trUseGMRedi) ) THEN
488           fMer(i,j) = fMer(i,j) + df(i,j)            diagName = 'DFyE'//diagSufx
489          ENDDO            CALL DIAGNOSTICS_FILL( df, diagName, k,1, 2,bi,bj, myThid )
        ENDDO  
490        ENDIF        ENDIF
491    #endif
 #ifdef NONLIN_FRSURF  
 C--   Compute vertical flux fVerT(kDown) at interface k+1 (between k & k+1):  
       IF ( calcAdvection .AND. K.EQ.Nr .AND.  
      &     useRealFreshWaterFlux .AND.  
      &     buoyancyRelation .EQ. 'OCEANICP' ) THEN    
        DO j=1-Oly,sNy+Oly  
         DO i=1-Olx,sNx+Olx  
          fVerT(i,j,kDown) = convertEmP2rUnit*PmEpR(i,j,bi,bj)  
      &     *rA(i,j,bi,bj)*maskC(i,j,k,bi,bj)*Tracer(i,j,k,bi,bj)  
         ENDDO  
        ENDDO  
       ENDIF  
 #endif /* NONLIN_FRSURF */  
492    
493  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):
494  C-    Advective flux in R  C-    Advective flux in R
495        IF (calcAdvection) THEN  #ifdef ALLOW_AIM
496  C     Note: wVel needs to be masked  C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr
497        IF (K.GE.2) THEN        IF (calcAdvection .AND. .NOT.implicitAdvection .AND. k.GE.2 .AND.
498         &     (.NOT.useAIM .OR. trIdentity.NE.GAD_SALINITY .OR. k.LT.Nr)
499         &   ) THEN
500    #else
501          IF (calcAdvection .AND. .NOT.implicitAdvection .AND. k.GE.2) THEN
502    #endif
503  C-    Compute vertical advective flux in the interior:  C-    Compute vertical advective flux in the interior:
504         IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN          IF (vertAdvecScheme.EQ.ENUM_CENTERED_2ND) THEN
505          CALL GAD_C2_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)             CALL GAD_C2_ADV_R( bi,bj,k, rTrans, TracAB, af, myThid )
506         ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF ( vertAdvecScheme.EQ.ENUM_UPWIND_1RST
507          CALL GAD_FLUXLIMIT_ADV_R(       &          .OR. vertAdvecScheme.EQ.ENUM_DST2 ) THEN
508       &       bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)             CALL GAD_DST2U1_ADV_R( bi,bj,k,vertAdvecScheme,deltaTLev(k),
509         ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN       I              rTrans, wFld, TracAB(1-OLx,1-OLy,1,bi,bj),
510          CALL GAD_U3_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)       O              af, myThid )
511         ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN          ELSEIF (vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN
512          CALL GAD_C4_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)             CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k, deltaTLev(k),
513         ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN       I              rTrans, wFld, TracAB(1-OLx,1-OLy,1,bi,bj),
514          CALL GAD_DST3_ADV_R(       O              af, myThid )
515       &       bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)          ELSEIF (vertAdvecScheme.EQ.ENUM_UPWIND_3RD ) THEN
516         ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN             CALL GAD_U3_ADV_R( bi,bj,k, rTrans, TracAB, af, myThid )
517          CALL GAD_DST3FL_ADV_R(          ELSEIF (vertAdvecScheme.EQ.ENUM_CENTERED_4TH) THEN
518       &       bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)             CALL GAD_C4_ADV_R( bi,bj,k, rTrans, TracAB, af, myThid )
519         ELSE          ELSEIF (vertAdvecScheme.EQ.ENUM_DST3 ) THEN
520          STOP 'GAD_CALC_RHS: Bad advectionScheme (R)'             CALL GAD_DST3_ADV_R( bi,bj,k, deltaTLev(k),
521         ENDIF       I              rTrans, wFld, TracAB(1-OLx,1-OLy,1,bi,bj),
522  C-    Surface "correction" term at k>1 :       O              af, myThid )
523         DO j=1-Oly,sNy+Oly          ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
524          DO i=1-Olx,sNx+Olx  cph This block is to trick the adjoint:
525           af(i,j) = af(i,j)  cph IF inAdExact=.FALSE., we want to use DST3
526       &           + (maskC(i,j,k,bi,bj)-maskC(i,j,k-1,bi,bj))*  cph with limiters in forward, but without limiters in reverse.
527       &             rTrans(i,j)*Tracer(i,j,k,bi,bj)            IF ( inAdMode ) THEN
528          ENDDO             CALL GAD_DST3_ADV_R( bi,bj,k, deltaTLev(k),
529         ENDDO       I              rTrans, wFld, TracAB(1-OLx,1-OLy,1,bi,bj),
530        ELSE       O              af, myThid )
531  C-    Surface "correction" term at k=1 :            ELSE
532         DO j=1-Oly,sNy+Oly             CALL GAD_DST3FL_ADV_R( bi,bj,k, deltaTLev(k),
533          DO i=1-Olx,sNx+Olx       I              rTrans, wFld, TracAB(1-OLx,1-OLy,1,bi,bj),
534           af(i,j) = rTrans(i,j)*Tracer(i,j,k,bi,bj)       O              af, myThid )
535          ENDDO            ENDIF
536         ENDDO  #ifndef ALLOW_AUTODIFF_TAMC
537        ENDIF          ELSEIF (vertAdvecScheme.EQ.ENUM_OS7MP ) THEN
538  C-    add the advective flux to fVerT             CALL GAD_OS7MP_ADV_R( bi,bj,k, deltaTLev(k),
539        DO j=1-Oly,sNy+Oly       I              rTrans, wFld, TracAB(1-OLx,1-OLy,1,bi,bj),
540         DO i=1-Olx,sNx+Olx       O              af, myThid )
541          fVerT(i,j,kUp) = fVerT(i,j,kUp) + af(i,j)  #endif
542         ENDDO          ELSE
543        ENDDO            STOP 'GAD_CALC_RHS: Bad vertAdvecScheme (R)'
544            ENDIF
545    C-     add the advective flux to fVerT
546            DO j=1-OLy,sNy+OLy
547             DO i=1-OLx,sNx+OLx
548              fVerT(i,j,kUp) = fVerT(i,j,kUp) + af(i,j)*maskInC(i,j,bi,bj)
549             ENDDO
550            ENDDO
551    #ifdef ALLOW_DIAGNOSTICS
552            IF ( useDiagnostics ) THEN
553              diagName = 'ADVr'//diagSufx
554              CALL DIAGNOSTICS_FILL( af, diagName, k,1, 2,bi,bj, myThid )
555    C- note: needs to explicitly increment the counter since DIAGNOSTICS_FILL
556    C        does it only if k=1 (never the case here)
557              IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT(diagName,bi,bj,myThid)
558            ENDIF
559    #endif
560        ENDIF        ENDIF
561    
562  C-    Diffusive flux in R  C-    Diffusive flux in R
563  C     Note: For K=1 then KM1=1 and this gives a dT/dr = 0 upper  C     Note: For K=1 then KM1=1 and this gives a dT/dr = 0 upper
564  C           boundary condition.  C           boundary condition.
565        IF (implicitDiffusion) THEN        IF (implicitDiffusion) THEN
566         DO j=1-Oly,sNy+Oly         DO j=1-OLy,sNy+OLy
567          DO i=1-Olx,sNx+Olx          DO i=1-OLx,sNx+OLx
568           df(i,j) = 0. _d 0           df(i,j) = 0. _d 0
569          ENDDO          ENDDO
570         ENDDO         ENDDO
571        ELSE        ELSE
572         CALL GAD_DIFF_R(bi,bj,k,KappaRT,tracer,df,myThid)         IF ( applyAB_onTracer ) THEN
573             CALL GAD_DIFF_R(bi,bj,k,KappaR,TracerN,df,myThid)
574           ELSE
575             CALL GAD_DIFF_R(bi,bj,k,KappaR,TracAB, df,myThid)
576           ENDIF
577          ENDIF
578    
579          IF ( trUseDiffKr4 ) THEN
580           IF ( applyAB_onTracer ) THEN
581             CALL GAD_BIHARM_R( bi,bj,k, diffKr4, TracerN, df, myThid )
582           ELSE
583             CALL GAD_BIHARM_R( bi,bj,k, diffKr4, TracAB,  df, myThid )
584           ENDIF
585        ENDIF        ENDIF
586    
587  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
588  C-    GM/Redi flux in R  C-    GM/Redi flux in R
589        IF (useGMRedi) THEN        IF ( trUseGMRedi ) THEN
590  C *note* should update GMREDI_RTRANSPORT to set df  *aja*  C *note* should update GMREDI_RTRANSPORT to set df  *aja*
591         CALL GMREDI_RTRANSPORT(          IF ( applyAB_onTracer ) THEN
592       I     iMin,iMax,jMin,jMax,bi,bj,K,            CALL GMREDI_RTRANSPORT(
593       I     Tracer,tracerIdentity,       I         iMin,iMax,jMin,jMax,bi,bj,k,
594       U     df,       I         TracerN,trIdentity,
595       I     myThid)       U         df,
596         I         myThid)
597            ELSE
598              CALL GMREDI_RTRANSPORT(
599         I         iMin,iMax,jMin,jMax,bi,bj,k,
600         I         TracAB, trIdentity,
601         U         df,
602         I         myThid)
603            ENDIF
604        ENDIF        ENDIF
605  #endif  #endif
606    
607        DO j=1-Oly,sNy+Oly        DO j=1-OLy,sNy+OLy
608         DO i=1-Olx,sNx+Olx         DO i=1-OLx,sNx+OLx
609          fVerT(i,j,kUp) = fVerT(i,j,kUp) + df(i,j)*maskUp(i,j)          fVerT(i,j,kUp) = fVerT(i,j,kUp) + df(i,j)*maskUp(i,j)
610         ENDDO         ENDDO
611        ENDDO        ENDDO
612    
613    #ifdef ALLOW_DIAGNOSTICS
614    C-    Diagnostics of Tracer flux in R dir (mainly Diffusive terms),
615    C       Explicit terms only & excluding advective terms:
616          IF ( useDiagnostics .AND.
617         &    (.NOT.implicitDiffusion .OR. trUseGMRedi) ) THEN
618              diagName = 'DFrE'//diagSufx
619              CALL DIAGNOSTICS_FILL( df, diagName, k,1, 2,bi,bj, myThid )
620          ENDIF
621    #endif
622    
623  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
624  C-    Add non local KPP transport term (ghat) to diffusive T flux.  C-    Set non local KPP transport term (ghat):
625        IF (useKPP) THEN        IF ( trUseKPP .AND. k.GE.2 ) THEN
626         DO j=1-Oly,sNy+Oly         DO j=1-OLy,sNy+OLy
627          DO i=1-Olx,sNx+Olx          DO i=1-OLx,sNx+OLx
628           df(i,j) = 0. _d 0           df(i,j) = 0. _d 0
629          ENDDO          ENDDO
630         ENDDO         ENDDO
631         IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN         IF (trIdentity.EQ.GAD_TEMPERATURE) THEN
 C *note* should update KPP_TRANSPORT_T to set df  *aja*  
632          CALL KPP_TRANSPORT_T(          CALL KPP_TRANSPORT_T(
633       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,       I           iMin,iMax,jMin,jMax,bi,bj,k,km1,
634       I     KappaRT,       O           df,
635       U     df )       I           myTime, myIter, myThid )
636         ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN         ELSEIF (trIdentity.EQ.GAD_SALINITY) THEN
637          CALL KPP_TRANSPORT_S(          CALL KPP_TRANSPORT_S(
638       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,       I           iMin,iMax,jMin,jMax,bi,bj,k,km1,
639       I     KappaRT,       O           df,
640       U     df )       I           myTime, myIter, myThid )
641    #ifdef ALLOW_PTRACERS
642           ELSEIF (trIdentity .GE. GAD_TR1) THEN
643            CALL KPP_TRANSPORT_PTR(
644         I           iMin,iMax,jMin,jMax,bi,bj,k,km1,
645         I           trIdentity-GAD_TR1+1,
646         O           df,
647         I           myTime, myIter, myThid )
648    #endif
649         ELSE         ELSE
650          STOP 'GAD_CALC_RHS: Ooops'          WRITE(errorMessageUnit,*)
651         &    'tracer identity =', trIdentity, ' is not valid => STOP'
652            STOP 'ABNORMAL END: S/R GAD_CALC_RHS: invalid tracer identity'
653         ENDIF         ENDIF
654         DO j=1-Oly,sNy+Oly         DO j=1-OLy,sNy+OLy
655          DO i=1-Olx,sNx+Olx          DO i=1-OLx,sNx+OLx
656           fVerT(i,j,kUp) = fVerT(i,j,kUp) + df(i,j)*maskUp(i,j)           fVerT(i,j,kUp) = fVerT(i,j,kUp)
657         &                  + df(i,j)*maskUp(i,j)*rhoFacF(k)
658            ENDDO
659           ENDDO
660    #ifdef ALLOW_DIAGNOSTICS
661    C-    Diagnostics of Non-Local Tracer (vertical) flux
662           IF ( useDiagnostics ) THEN
663             diagName = 'KPPg'//diagSufx
664             CALL DIAGNOSTICS_FILL( df, diagName, k,1, 2,bi,bj, myThid )
665    C- note: needs to explicitly increment the counter since DIAGNOSTICS_FILL
666    C        does it only if k=1 (never the case here)
667             IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT(diagName,bi,bj,myThid)
668           ENDIF
669    #endif
670          ENDIF
671    #endif /* ALLOW_KPP */
672    
673    #ifdef GAD_SMOLARKIEWICZ_HACK
674    coj   Hack to make redi (and everything else in this s/r) positive
675    coj   (see Smolarkiewicz MWR 1989 and Bott MWR 1989).
676    coj   Only works if 'down' is k+1 and k loop in thermodynamics is k=Nr,1,-1
677    coj
678    coj   Apply to all tracers except temperature
679          IF ( trIdentity.NE.GAD_TEMPERATURE .AND.
680         &     trIdentity.NE.GAD_SALINITY ) THEN
681           DO j=1-OLy,sNy+OLy-1
682            DO i=1-OLx,sNx+OLx-1
683    coj   Add outgoing fluxes
684             outFlux=deltaTLev(k)*
685         &    _recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
686         &   *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)*recip_rhoFacC(k)
687         &    *( MAX(0. _d 0,fZon(i+1,j)) + MAX(0. _d 0,-fZon(i,j))
688         &      +MAX(0. _d 0,fMer(i,j+1)) + MAX(0. _d 0,-fMer(i,j))
689         &      +MAX(0. _d 0,fVerT(i,j,kDown)*rkSign)
690         &      +MAX(0. _d 0,-fVerT(i,j,kUp)*rkSign)
691         &     )
692             IF ( applyAB_onTracer ) THEN
693               trac=TracerN(i,j,k,bi,bj)
694             ELSE
695               trac=TracAB(i,j,k,bi,bj)
696             ENDIF
697    coj   If they would reduce tracer by a fraction of more than
698    coj   SmolarkiewiczMaxFrac, scale them down
699             IF (outFlux.GT.0. _d 0 .AND.
700         &       outFlux.GT.SmolarkiewiczMaxFrac*trac) THEN
701    coj   If tracer is already negative, scale flux to zero
702               fac = MAX(0. _d 0,SmolarkiewiczMaxFrac*trac/outFlux)
703    
704               IF (fZon(i+1,j).GT.0. _d 0) fZon(i+1,j)=fac*fZon(i+1,j)
705               IF (-fZon(i,j) .GT.0. _d 0) fZon(i,j)  =fac*fZon(i,j)
706               IF (fMer(i,j+1).GT.0. _d 0) fMer(i,j+1)=fac*fMer(i,j+1)
707               IF (-fMer(i,j) .GT.0. _d 0) fMer(i,j)  =fac*fMer(i,j)
708               IF (-fVerT(i,j,kUp)*rkSign .GT.0. _d 0)
709         &       fVerT(i,j,kUp)=fac*fVerT(i,j,kUp)
710    
711               IF (k.LT.Nr .AND. fVerT(i,j,kDown)*rkSign.GT.0. _d 0) THEN
712    coj   Down flux is special: it has already been applied in lower layer,
713    coj   so we have to readjust this.
714    coj   Note: for k+1, gTracer is now the updated tracer, not the tendency!
715    coj   thus it has an extra factor deltaTLev(k+1)
716                 gTrFac=deltaTLev(k+1)
717    coj   Other factors that have been applied to gTracer since the last call:
718    #ifdef NONLIN_FRSURF
719                 IF (nonlinFreeSurf.GT.0) THEN
720                  IF (select_rStar.GT.0) THEN
721    #ifndef DISABLE_RSTAR_CODE
722                    gTrFac = gTrFac/rStarExpC(i,j,bi,bj)
723    #endif /* DISABLE_RSTAR_CODE */
724                  ENDIF
725                 ENDIF
726    #endif /* NONLIN_FRSURF */
727    coj   Now: undo down flux, ...
728                 gTracer(i,j,k+1,bi,bj)=gTracer(i,j,k+1,bi,bj)
729         &        +gTrFac
730         &         *_recip_hFacC(i,j,k+1,bi,bj)*recip_drF(k+1)
731         &         *recip_rA(i,j,bi,bj)*recip_deepFac2C(k+1)
732         &         *recip_rhoFacC(k+1)
733         &         *( -fVerT(i,j,kDown)*rkSign )
734    coj   ... scale ...
735                 fVerT(i,j,kDown)=fac*fVerT(i,j,kDown)
736    coj   ... and reapply
737                 gTracer(i,j,k+1,bi,bj)=gTracer(i,j,k+1,bi,bj)
738         &        +gTrFac
739         &         *_recip_hFacC(i,j,k+1,bi,bj)*recip_drF(k+1)
740         &         *recip_rA(i,j,bi,bj)*recip_deepFac2C(k+1)
741         &         *recip_rhoFacC(k+1)
742         &         *( fVerT(i,j,kDown)*rkSign )
743               ENDIF
744    
745             ENDIF
746          ENDDO          ENDDO
747         ENDDO         ENDDO
748        ENDIF        ENDIF
749  #endif  #endif
750    
751  C--   Divergence of fluxes  C--   Divergence of fluxes
752        DO j=1-Oly,sNy+Oly-1  C     Anelastic: scale vertical fluxes by rhoFac and leave Horizontal fluxes unchanged
753         DO i=1-Olx,sNx+Olx-1  C     for Stevens OBC: keep only vertical diffusive contribution on boundaries
754          DO j=1-OLy,sNy+OLy-1
755           DO i=1-OLx,sNx+OLx-1
756          gTracer(i,j,k,bi,bj)=gTracer(i,j,k,bi,bj)          gTracer(i,j,k,bi,bj)=gTracer(i,j,k,bi,bj)
757       &   -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &   -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
758       &    *recip_rA(i,j,bi,bj)       &   *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)*recip_rhoFacC(k)
759       &    *(       &   *( (fZon(i+1,j)-fZon(i,j))*maskInC(i,j,bi,bj)
760       &    +( fZon(i+1,j)-fZon(i,j) )       &     +(fMer(i,j+1)-fMer(i,j))*maskInC(i,j,bi,bj)
761       &    +( fMer(i,j+1)-fMer(i,j) )       &     +(fVerT(i,j,kDown)-fVerT(i,j,kUp))*rkSign
762       &    +( fVerT(i,j,kUp)-fVerT(i,j,kDown) )*rkFac       &     -localT(i,j)*( (uTrans(i+1,j)-uTrans(i,j))*advFac
763         &                   +(vTrans(i,j+1)-vTrans(i,j))*advFac
764         &                   +(rTransKp1(i,j)-rTrans(i,j))*rAdvFac
765         &                  )*maskInC(i,j,bi,bj)
766       &    )       &    )
767         ENDDO         ENDDO
768        ENDDO        ENDDO
769    
770  #ifdef NONLIN_FRSURF  #ifdef ALLOW_DEBUG
771  C-- account for 3.D divergence of the flow in rStar coordinate:        IF ( debugLevel .GE. debLevC
772        IF (calcAdvection .AND. select_rStar.GT.0) THEN       &   .AND. trIdentity.EQ.GAD_TEMPERATURE
773         DO j=1-Oly,sNy+Oly-1       &   .AND. k.EQ.2 .AND. myIter.EQ.1+nIter0
774          DO i=1-Olx,sNx+Olx-1       &   .AND. nPx.EQ.1 .AND. nPy.EQ.1
775           gTracer(i,j,k,bi,bj) = gTracer(i,j,k,bi,bj)       &   .AND. useCubedSphereExchange ) THEN
776       &     - (rStarExpC(i,j,bi,bj) - 1. _d 0)/deltaTfreesurf          CALL DEBUG_CS_CORNER_UV( ' fZon,fMer from GAD_CALC_RHS',
777       &       *tracer(i,j,k,bi,bj)*maskC(i,j,k,bi,bj)       &             fZon,fMer, k, standardMessageUnit,bi,bj,myThid )
         ENDDO  
        ENDDO  
       ENDIF  
       IF (calcAdvection .AND. select_rStar.LT.0) THEN  
        DO j=1-Oly,sNy+Oly-1  
         DO i=1-Olx,sNx+Olx-1  
          gTracer(i,j,k,bi,bj) = gTracer(i,j,k,bi,bj)  
      &     - rStarDhCDt(i,j,bi,bj)  
      &       *tracer(i,j,k,bi,bj)*maskC(i,j,k,bi,bj)  
         ENDDO  
        ENDDO  
778        ENDIF        ENDIF
779  #endif /* NONLIN_FRSURF */  #endif /* ALLOW_DEBUG */
         
780    
781        RETURN        RETURN
782        END        END
Legend:
Removed from v.1.17  
changed lines
  Added in v.1.61
  ViewVC Help
Powered by ViewVC 1.1.22