/[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.42 by jmc, Sun Nov 19 22:04:09 2006 UTC revision 1.61 by jmc, Fri Jun 15 13:28:41 2012 UTC
# Line 11  C !INTERFACE: ========================== Line 11  C !INTERFACE: ==========================
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, maskUp, uFld, vFld, wFld,       I           xA, yA, maskUp, uFld, vFld, wFld,
13       I           uTrans, vTrans, rTrans, rTransKp1,       I           uTrans, vTrans, rTrans, rTransKp1,
14       I           diffKh, diffK4, KappaR, TracerN, TracAB,       I           diffKh, diffK4, KappaR, diffKr4, TracerN, TracAB,
15       I           tracerIdentity, advectionScheme, vertAdvecScheme,       I           deltaTLev, trIdentity,
16         I           advectionScheme, vertAdvecScheme,
17       I           calcAdvection, implicitAdvection, applyAB_onTracer,       I           calcAdvection, implicitAdvection, applyAB_onTracer,
18         I           trUseDiffKr4, trUseGMRedi, trUseKPP,
19       U           fVerT, gTracer,       U           fVerT, gTracer,
20       I           myTime, myIter, myThid )       I           myTime, myIter, myThid )
21    
# Line 65  C uTrans,vTrans    :: 2-D arrays of volu Line 67  C uTrans,vTrans    :: 2-D arrays of volu
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
69  C diffKh           :: horizontal diffusion coefficient  C diffKh           :: horizontal diffusion coefficient
70  C diffK4           :: bi-harmonic diffusion coefficient  C diffK4           :: horizontal 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 diffKr4          :: 1-D array for vertical bi-harmonic diffusion coefficient
73  C TracerN          :: tracer field @ time-step n (Note: only used  C TracerN          :: tracer field @ time-step n (Note: only used
74  C                     if applying AB on tracer field rather than on tendency gTr)  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  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)  C                     or extrapolated fwd in time to n+1/2 if applying AB on Tr)
77  C tracerIdentity   :: tracer identifier (required for KPP,GM)  C trIdentity       :: tracer identifier (required for KPP,GM)
78  C advectionScheme  :: advection scheme to use (Horizontal plane)  C advectionScheme  :: advection scheme to use (Horizontal plane)
79  C vertAdvecScheme  :: advection scheme to use (Vertical direction)  C vertAdvecScheme  :: advection scheme to use (Vertical direction)
80  C calcAdvection    :: =False if Advec computed with multiDim scheme  C calcAdvection    :: =False if Advec computed with multiDim scheme
81  C implicitAdvection:: =True if vertical Advec computed implicitly  C implicitAdvection:: =True if vertical Advec computed implicitly
82  C applyAB_onTracer :: apply Adams-Bashforth on Tracer (rather than on gTr)  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  C myTime           :: current time
87  C myIter           :: iteration number  C myIter           :: iteration number
88  C myThid           :: thread number  C myThid           :: thread number
# Line 94  C myThid           :: thread number Line 100  C myThid           :: thread number
100        _RL rTransKp1(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 KappaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL KappaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
103          _RL diffKr4(Nr)
104        _RL TracerN(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL TracerN(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
105        _RL TracAB (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL TracAB (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
106        INTEGER tracerIdentity        _RL deltaTLev(Nr)
107          INTEGER trIdentity
108        INTEGER advectionScheme, vertAdvecScheme        INTEGER advectionScheme, vertAdvecScheme
109        LOGICAL calcAdvection        LOGICAL calcAdvection
110        LOGICAL implicitAdvection, applyAB_onTracer        LOGICAL implicitAdvection, applyAB_onTracer
111          LOGICAL trUseDiffKr4, trUseGMRedi, trUseKPP
112        _RL     myTime        _RL     myTime
113        INTEGER myIter, myThid        INTEGER myIter, myThid
114    
# Line 109  C fVerT            :: 2 1/2D arrays for Line 118  C fVerT            :: 2 1/2D arrays for
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
# Line 118  C af               :: advective flux Line 133  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  C locABT           :: local copy of (AB-extrapolated) tracer field
 #ifdef ALLOW_DIAGNOSTICS  
       CHARACTER*8 diagName  
       CHARACTER*4 GAD_DIAG_SUFX, diagSufx  
       EXTERNAL    GAD_DIAG_SUFX  
 #endif  
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)
# Line 132  C locABT           :: local copy of (AB- Line 144  C locABT           :: local copy of (AB-
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)        _RL locABT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
146        _RL advFac, rAdvFac        _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 143  C--   the kDown is still required Line 162  C--   the kDown is still required
162  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
163  C--   Set diagnostic suffix for the current tracer  C--   Set diagnostic suffix for the current tracer
164        IF ( useDiagnostics ) THEN        IF ( useDiagnostics ) THEN
165          diagSufx = GAD_DIAG_SUFX( tracerIdentity, myThid )          diagSufx = GAD_DIAG_SUFX( trIdentity, myThid )
166        ENDIF        ENDIF
167  #endif  #endif
168    
169        advFac  = 0. _d 0        advFac  = 0. _d 0
170        IF (calcAdvection) advFac = 1. _d 0        IF (calcAdvection) advFac = 1. _d 0
171        rAdvFac = rkSign*advFac        rAdvFac = rkSign*advFac
172        IF (implicitAdvection) rAdvFac = 0. _d 0        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
# Line 183  C--   Unless we have already calculated Line 202  C--   Unless we have already calculated
202  C     the tendency to zero.  C     the tendency to zero.
203  C     <== now done earlier at the beginning of thermodynamics.  C     <== now done earlier at the beginning of thermodynamics.
204  c     IF (calcAdvection) THEN  c     IF (calcAdvection) THEN
205  c      DO j=1-Oly,sNy+Oly  c      DO j=1-OLy,sNy+OLy
206  c       DO i=1-Olx,sNx+Olx  c       DO i=1-OLx,sNx+OLx
207  c        gTracer(i,j,k,bi,bj)=0. _d 0  c        gTracer(i,j,k,bi,bj)=0. _d 0
208  c       ENDDO  c       ENDDO
209  c      ENDDO  c      ENDDO
# Line 198  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
# Line 207  C--   Initialize net flux in X direction Line 226  C--   Initialize net flux in X direction
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,locABT,af,myThid)             CALL GAD_C2_ADV_X( bi,bj,k, uTrans, locABT, af, myThid )
230          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST
231       &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN       &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
232            CALL GAD_DST2U1_ADV_X( bi,bj,k, advectionScheme,             CALL GAD_DST2U1_ADV_X( bi,bj,k, advectionScheme, .TRUE.,
233       I            dTtracerLev(k), uTrans, uFld, locABT,       I              deltaTLev(k), uTrans, uFld, locABT,
234       O            af, myThid )       O              af, myThid )
235          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSE
236            CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, dTtracerLev(k),           DO j=1-OLy,sNy+OLy
237       I            uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,            DO i=1-OLx,sNx+OLx
238       O            af, myThid )  #ifdef ALLOW_OBCS
239          ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN             maskLocW(i,j) = _maskW(i,j,k,bi,bj)*maskInW(i,j,bi,bj)
240            CALL GAD_U3_ADV_X(bi,bj,k,uTrans,locABT,af,myThid)  #else /* ALLOW_OBCS */
241          ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN             maskLocW(i,j) = _maskW(i,j,k,bi,bj)
242            CALL GAD_C4_ADV_X(bi,bj,k,uTrans,locABT,af,myThid)  #endif /* ALLOW_OBCS */
243          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN            ENDDO
244            CALL GAD_DST3_ADV_X( bi,bj,k, dTtracerLev(k),           ENDDO
245       I            uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,           IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
246       O            af, myThid )             CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
247          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN       I              uTrans, uFld, maskLocW, locABT,
248           IF ( inAdMode ) THEN       O              af, myThid )
249             ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
250               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:  cph This block is to trick the adjoint:
262  cph IF inAdExact=.FALSE., we want to use DST3  cph IF inAdExact=.FALSE., we want to use DST3
263  cph with limiters in forward, but without limiters in reverse.  cph with limiters in forward, but without limiters in reverse.
264            CALL GAD_DST3_ADV_X( bi,bj,k, dTtracerLev(k),             CALL GAD_DST3_ADV_X( bi,bj,k, .TRUE., deltaTLev(k),
265       I           uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,       I              uTrans, uFld, maskLocW, locABT,
266       O           af, myThid )       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           ELSE
279            CALL GAD_DST3FL_ADV_X( bi,bj,k, dTtracerLev(k),            STOP 'GAD_CALC_RHS: Bad advectionScheme (X)'
      I           uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,  
      O           af, myThid )  
280           ENDIF           ENDIF
         ELSE  
          STOP 'GAD_CALC_RHS: Bad advectionScheme (X)'  
281          ENDIF          ENDIF
282          DO j=1-Oly,sNy+Oly  #ifdef ALLOW_OBCS
283           DO i=1-Olx,sNx+Olx          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)            fZon(i,j) = fZon(i,j) + af(i,j)
294           ENDDO           ENDDO
295          ENDDO          ENDDO
296  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
297          IF ( useDiagnostics ) THEN          IF ( useDiagnostics ) THEN
298            diagName = 'ADVx'//diagSufx            diagName = 'ADVx'//diagSufx
299            CALL DIAGNOSTICS_FILL(af,diagName, k,1, 2,bi,bj, myThid)            CALL DIAGNOSTICS_FILL( af, diagName, k,1, 2,bi,bj, myThid )
300          ENDIF          ENDIF
301  #endif  #endif
302        ENDIF        ENDIF
# Line 258  C-    Diffusive flux in X Line 305  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
# Line 272  C-    Add bi-harmonic diffusive flux in Line 319  C-    Add bi-harmonic diffusive flux in
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 set df  *aja*  C *note* should update GMREDI_XTRANSPORT to set df  *aja*
324          IF ( applyAB_onTracer ) THEN          IF ( applyAB_onTracer ) THEN
325            CALL GMREDI_XTRANSPORT(            CALL GMREDI_XTRANSPORT(
326       I         iMin,iMax,jMin,jMax,bi,bj,k,       I         iMin,iMax,jMin,jMax,bi,bj,k,
327       I         xA,TracerN,tracerIdentity,       I         xA,TracerN,trIdentity,
328       U         df,       U         df,
329       I         myThid)       I         myThid)
330          ELSE          ELSE
331            CALL GMREDI_XTRANSPORT(            CALL GMREDI_XTRANSPORT(
332       I         iMin,iMax,jMin,jMax,bi,bj,k,       I         iMin,iMax,jMin,jMax,bi,bj,k,
333       I         xA,TracAB, tracerIdentity,       I         xA,TracAB, trIdentity,
334       U         df,       U         df,
335       I         myThid)       I         myThid)
336          ENDIF          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    
# Line 299  C *note* should update GMREDI_XTRANSPORT Line 347  C *note* should update GMREDI_XTRANSPORT
347  C-    Diagnostics of Tracer flux in X dir (mainly Diffusive term),  C-    Diagnostics of Tracer flux in X dir (mainly Diffusive term),
348  C       excluding advective terms:  C       excluding advective terms:
349        IF ( useDiagnostics .AND.        IF ( useDiagnostics .AND.
350       &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. useGMRedi) ) THEN       &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. trUseGMRedi) ) THEN
351            diagName = 'DFxE'//diagSufx            diagName = 'DFxE'//diagSufx
352            CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)            CALL DIAGNOSTICS_FILL( df, diagName, k,1, 2,bi,bj, myThid )
353        ENDIF        ENDIF
354  #endif  #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
# Line 315  C--   Initialize net flux in Y direction Line 363  C--   Initialize net flux in Y direction
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,locABT,af,myThid)             CALL GAD_C2_ADV_Y( bi,bj,k, vTrans, locABT, af, myThid )
367          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST
368       &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN       &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
369            CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme,             CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme, .TRUE.,
370       I            dTtracerLev(k), vTrans, vFld, locABT,       I              deltaTLev(k), vTrans, vFld, locABT,
371       O            af, myThid )       O              af, myThid )
372          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSE
373            CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, dTtracerLev(k),           DO j=1-OLy,sNy+OLy
374       I            vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,            DO i=1-OLx,sNx+OLx
375       O            af, myThid )  #ifdef ALLOW_OBCS
376          ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN             maskLocS(i,j) = _maskS(i,j,k,bi,bj)*maskInS(i,j,bi,bj)
377            CALL GAD_U3_ADV_Y(bi,bj,k,vTrans,locABT,af,myThid)  #else /* ALLOW_OBCS */
378          ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN             maskLocS(i,j) = _maskS(i,j,k,bi,bj)
379            CALL GAD_C4_ADV_Y(bi,bj,k,vTrans,locABT,af,myThid)  #endif /* ALLOW_OBCS */
380          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN            ENDDO
381            CALL GAD_DST3_ADV_Y( bi,bj,k, dTtracerLev(k),           ENDDO
382       I            vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,           IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
383       O            af, myThid )             CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
384          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN       I              vTrans, vFld, maskLocS, locABT,
385           IF ( inAdMode ) THEN       O              af, myThid )
386             ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
387               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:  cph This block is to trick the adjoint:
399  cph IF inAdExact=.FALSE., we want to use DST3  cph IF inAdExact=.FALSE., we want to use DST3
400  cph with limiters in forward, but without limiters in reverse.  cph with limiters in forward, but without limiters in reverse.
401            CALL GAD_DST3_ADV_Y( bi,bj,k, dTtracerLev(k),             CALL GAD_DST3_ADV_Y( bi,bj,k, .TRUE., deltaTLev(k),
402       I           vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,       I              vTrans, vFld, maskLocS, locABT,
403       O           af, myThid )       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           ELSE
416            CALL GAD_DST3FL_ADV_Y( bi,bj,k, dTtracerLev(k),             STOP 'GAD_CALC_RHS: Bad advectionScheme (Y)'
      I           vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,  
      O           af, myThid )  
417           ENDIF           ENDIF
         ELSE  
           STOP 'GAD_CALC_RHS: Bad advectionScheme (Y)'  
418          ENDIF          ENDIF
419          DO j=1-Oly,sNy+Oly  #ifdef ALLOW_OBCS
420           DO i=1-Olx,sNx+Olx          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)            fMer(i,j) = fMer(i,j) + af(i,j)
431           ENDDO           ENDDO
432          ENDDO          ENDDO
433  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
434          IF ( useDiagnostics ) THEN          IF ( useDiagnostics ) THEN
435            diagName = 'ADVy'//diagSufx            diagName = 'ADVy'//diagSufx
436            CALL DIAGNOSTICS_FILL(af,diagName, k,1, 2,bi,bj, myThid)            CALL DIAGNOSTICS_FILL( af, diagName, k,1, 2,bi,bj, myThid )
437          ENDIF          ENDIF
438  #endif  #endif
439        ENDIF        ENDIF
# Line 366  C-    Diffusive flux in Y Line 442  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
# Line 380  C-    Add bi-harmonic flux in Y Line 456  C-    Add bi-harmonic flux in Y
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 set df  *aja*  C *note* should update GMREDI_YTRANSPORT to set df  *aja*
461          IF ( applyAB_onTracer ) THEN          IF ( applyAB_onTracer ) THEN
462            CALL GMREDI_YTRANSPORT(            CALL GMREDI_YTRANSPORT(
463       I         iMin,iMax,jMin,jMax,bi,bj,k,       I         iMin,iMax,jMin,jMax,bi,bj,k,
464       I         yA,TracerN,tracerIdentity,       I         yA,TracerN,trIdentity,
465       U         df,       U         df,
466       I         myThid)       I         myThid)
467          ELSE          ELSE
468            CALL GMREDI_YTRANSPORT(            CALL GMREDI_YTRANSPORT(
469       I         iMin,iMax,jMin,jMax,bi,bj,k,       I         iMin,iMax,jMin,jMax,bi,bj,k,
470       I         yA,TracAB, tracerIdentity,       I         yA,TracAB, trIdentity,
471       U         df,       U         df,
472       I         myThid)       I         myThid)
473          ENDIF          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    
# Line 407  C *note* should update GMREDI_YTRANSPORT Line 484  C *note* should update GMREDI_YTRANSPORT
484  C-    Diagnostics of Tracer flux in Y dir (mainly Diffusive terms),  C-    Diagnostics of Tracer flux in Y dir (mainly Diffusive terms),
485  C       excluding advective terms:  C       excluding advective terms:
486        IF ( useDiagnostics .AND.        IF ( useDiagnostics .AND.
487       &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. useGMRedi) ) THEN       &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. trUseGMRedi) ) THEN
488            diagName = 'DFyE'//diagSufx            diagName = 'DFyE'//diagSufx
489            CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)            CALL DIAGNOSTICS_FILL( df, diagName, k,1, 2,bi,bj, myThid )
490        ENDIF        ENDIF
491  #endif  #endif
492    
# Line 418  C-    Advective flux in R Line 495  C-    Advective flux in R
495  #ifdef ALLOW_AIM  #ifdef ALLOW_AIM
496  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
497        IF (calcAdvection .AND. .NOT.implicitAdvection .AND. k.GE.2 .AND.        IF (calcAdvection .AND. .NOT.implicitAdvection .AND. k.GE.2 .AND.
498       &     (.NOT.useAIM .OR.tracerIdentity.NE.GAD_SALINITY .OR.k.LT.Nr)       &     (.NOT.useAIM .OR. trIdentity.NE.GAD_SALINITY .OR. k.LT.Nr)
499       &   ) THEN       &   ) THEN
500  #else  #else
501        IF (calcAdvection .AND. .NOT.implicitAdvection .AND. k.GE.2) THEN        IF (calcAdvection .AND. .NOT.implicitAdvection .AND. k.GE.2) THEN
502  #endif  #endif
503  C-    Compute vertical advective flux in the interior:  C-    Compute vertical advective flux in the interior:
504          IF (vertAdvecScheme.EQ.ENUM_CENTERED_2ND) THEN          IF (vertAdvecScheme.EQ.ENUM_CENTERED_2ND) THEN
505            CALL GAD_C2_ADV_R(bi,bj,k,rTrans,TracAB,af,myThid)             CALL GAD_C2_ADV_R( bi,bj,k, rTrans, TracAB, af, myThid )
506          ELSEIF ( vertAdvecScheme.EQ.ENUM_UPWIND_1RST          ELSEIF ( vertAdvecScheme.EQ.ENUM_UPWIND_1RST
507       &          .OR. vertAdvecScheme.EQ.ENUM_DST2 ) THEN       &          .OR. vertAdvecScheme.EQ.ENUM_DST2 ) THEN
508            CALL GAD_DST2U1_ADV_R( bi,bj,k, vertAdvecScheme,             CALL GAD_DST2U1_ADV_R( bi,bj,k,vertAdvecScheme,deltaTLev(k),
509       I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),       I              rTrans, wFld, TracAB(1-OLx,1-OLy,1,bi,bj),
510       O         af, myThid )       O              af, myThid )
511          ELSEIF (vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF (vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN
512            CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k,             CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k, deltaTLev(k),
513       I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),       I              rTrans, wFld, TracAB(1-OLx,1-OLy,1,bi,bj),
514       O         af, myThid )       O              af, myThid )
515          ELSEIF (vertAdvecScheme.EQ.ENUM_UPWIND_3RD ) THEN          ELSEIF (vertAdvecScheme.EQ.ENUM_UPWIND_3RD ) THEN
516            CALL GAD_U3_ADV_R(bi,bj,k,rTrans,TracAB,af,myThid)             CALL GAD_U3_ADV_R( bi,bj,k, rTrans, TracAB, af, myThid )
517          ELSEIF (vertAdvecScheme.EQ.ENUM_CENTERED_4TH) THEN          ELSEIF (vertAdvecScheme.EQ.ENUM_CENTERED_4TH) THEN
518            CALL GAD_C4_ADV_R(bi,bj,k,rTrans,TracAB,af,myThid)             CALL GAD_C4_ADV_R( bi,bj,k, rTrans, TracAB, af, myThid )
519          ELSEIF (vertAdvecScheme.EQ.ENUM_DST3 ) THEN          ELSEIF (vertAdvecScheme.EQ.ENUM_DST3 ) THEN
520            CALL GAD_DST3_ADV_R( bi,bj,k,             CALL GAD_DST3_ADV_R( bi,bj,k, deltaTLev(k),
521       I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),       I              rTrans, wFld, TracAB(1-OLx,1-OLy,1,bi,bj),
522       O         af, myThid )       O              af, myThid )
523          ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN          ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
524  cph This block is to trick the adjoint:  cph This block is to trick the adjoint:
525  cph IF inAdExact=.FALSE., we want to use DST3  cph IF inAdExact=.FALSE., we want to use DST3
526  cph with limiters in forward, but without limiters in reverse.  cph with limiters in forward, but without limiters in reverse.
527            IF ( inAdMode ) THEN            IF ( inAdMode ) THEN
528             CALL GAD_DST3_ADV_R( bi,bj,k,             CALL GAD_DST3_ADV_R( bi,bj,k, deltaTLev(k),
529       I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),       I              rTrans, wFld, TracAB(1-OLx,1-OLy,1,bi,bj),
530       O         af, myThid )       O              af, myThid )
531            ELSE            ELSE
532             CALL GAD_DST3FL_ADV_R( bi,bj,k,             CALL GAD_DST3FL_ADV_R( bi,bj,k, deltaTLev(k),
533       I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),       I              rTrans, wFld, TracAB(1-OLx,1-OLy,1,bi,bj),
534       O         af, myThid )       O              af, myThid )
535            ENDIF            ENDIF
536    #ifndef ALLOW_AUTODIFF_TAMC
537            ELSEIF (vertAdvecScheme.EQ.ENUM_OS7MP ) THEN
538               CALL GAD_OS7MP_ADV_R( bi,bj,k, deltaTLev(k),
539         I              rTrans, wFld, TracAB(1-OLx,1-OLy,1,bi,bj),
540         O              af, myThid )
541    #endif
542          ELSE          ELSE
543            STOP 'GAD_CALC_RHS: Bad vertAdvecScheme (R)'            STOP 'GAD_CALC_RHS: Bad vertAdvecScheme (R)'
544          ENDIF          ENDIF
545  C-     add the advective flux to fVerT  C-     add the advective flux to fVerT
546          DO j=1-Oly,sNy+Oly          DO j=1-OLy,sNy+OLy
547           DO i=1-Olx,sNx+Olx           DO i=1-OLx,sNx+OLx
548            fVerT(i,j,kUp) = fVerT(i,j,kUp) + af(i,j)            fVerT(i,j,kUp) = fVerT(i,j,kUp) + af(i,j)*maskInC(i,j,bi,bj)
549           ENDDO           ENDDO
550          ENDDO          ENDDO
551  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
552          IF ( useDiagnostics ) THEN          IF ( useDiagnostics ) THEN
553            diagName = 'ADVr'//diagSufx            diagName = 'ADVr'//diagSufx
554            CALL DIAGNOSTICS_FILL(af,diagName, k,1, 2,bi,bj, myThid)            CALL DIAGNOSTICS_FILL( af, diagName, k,1, 2,bi,bj, myThid )
555  C- note: needs to explicitly increment the counter since DIAGNOSTICS_FILL  C- note: needs to explicitly increment the counter since DIAGNOSTICS_FILL
556  C        does it only if k=1 (never the case here)  C        does it only if k=1 (never the case here)
557            IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT(diagName,bi,bj,myThid)            IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT(diagName,bi,bj,myThid)
# Line 480  C-    Diffusive flux in R Line 563  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
# Line 493  C           boundary condition. Line 576  C           boundary condition.
576         ENDIF         ENDIF
577        ENDIF        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
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          IF ( applyAB_onTracer ) THEN          IF ( applyAB_onTracer ) THEN
592            CALL GMREDI_RTRANSPORT(            CALL GMREDI_RTRANSPORT(
593       I         iMin,iMax,jMin,jMax,bi,bj,k,       I         iMin,iMax,jMin,jMax,bi,bj,k,
594       I         TracerN,tracerIdentity,       I         TracerN,trIdentity,
595       U         df,       U         df,
596       I         myThid)       I         myThid)
597          ELSE          ELSE
598            CALL GMREDI_RTRANSPORT(            CALL GMREDI_RTRANSPORT(
599       I         iMin,iMax,jMin,jMax,bi,bj,k,       I         iMin,iMax,jMin,jMax,bi,bj,k,
600       I         TracAB, tracerIdentity,       I         TracAB, trIdentity,
601       U         df,       U         df,
602       I         myThid)       I         myThid)
603          ENDIF          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
# Line 523  C *note* should update GMREDI_RTRANSPORT Line 614  C *note* should update GMREDI_RTRANSPORT
614  C-    Diagnostics of Tracer flux in R dir (mainly Diffusive terms),  C-    Diagnostics of Tracer flux in R dir (mainly Diffusive terms),
615  C       Explicit terms only & excluding advective terms:  C       Explicit terms only & excluding advective terms:
616        IF ( useDiagnostics .AND.        IF ( useDiagnostics .AND.
617       &    (.NOT.implicitDiffusion .OR. useGMRedi) ) THEN       &    (.NOT.implicitDiffusion .OR. trUseGMRedi) ) THEN
618            diagName = 'DFrE'//diagSufx            diagName = 'DFrE'//diagSufx
619            CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)            CALL DIAGNOSTICS_FILL( df, diagName, k,1, 2,bi,bj, myThid )
620        ENDIF        ENDIF
621  #endif  #endif
622    
623  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
624  C-    Set non local KPP transport term (ghat):  C-    Set non local KPP transport term (ghat):
625        IF ( useKPP .AND. k.GE.2 ) 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
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       O     df )       O           df,
635         ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN       I           myTime, myIter, myThid )
636           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       O     df )       O           df,
640         I           myTime, myIter, myThid )
641  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
642         ELSEIF (tracerIdentity .GE. GAD_TR1) THEN         ELSEIF (trIdentity .GE. GAD_TR1) THEN
643          CALL KPP_TRANSPORT_PTR(          CALL KPP_TRANSPORT_PTR(
644       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,       I           iMin,iMax,jMin,jMax,bi,bj,k,km1,
645       I     tracerIdentity-GAD_TR1+1,       I           trIdentity-GAD_TR1+1,
646       O     df )       O           df,
647         I           myTime, myIter, myThid )
648  #endif  #endif
649         ELSE         ELSE
650          PRINT*,'invalid tracer indentity: ', tracerIdentity          WRITE(errorMessageUnit,*)
651          STOP 'GAD_CALC_RHS: Ooops'       &    'tracer identity =', trIdentity, ' is not valid => STOP'
652            STOP 'ABNORMAL END: S/R GAD_CALC_RHS: invalid tracer identity'
653           ENDIF
654           DO j=1-OLy,sNy+OLy
655            DO i=1-OLx,sNx+OLx
656             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         ENDIF
669         DO j=1-Oly,sNy+Oly  #endif
670          DO i=1-Olx,sNx+Olx        ENDIF
671           fVerT(i,j,kUp) = fVerT(i,j,kUp) + df(i,j)*maskUp(i,j)  #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_rA(i,j,bi,bj)       &   -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
758       &   *( (fZon(i+1,j)-fZon(i,j))       &   *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)*recip_rhoFacC(k)
759       &     +(fMer(i,j+1)-fMer(i,j))       &   *( (fZon(i+1,j)-fZon(i,j))*maskInC(i,j,bi,bj)
760         &     +(fMer(i,j+1)-fMer(i,j))*maskInC(i,j,bi,bj)
761       &     +(fVerT(i,j,kDown)-fVerT(i,j,kUp))*rkSign       &     +(fVerT(i,j,kDown)-fVerT(i,j,kUp))*rkSign
762       &     -localT(i,j)*( (uTrans(i+1,j)-uTrans(i,j))       &     -localT(i,j)*( (uTrans(i+1,j)-uTrans(i,j))*advFac
763       &                   +(vTrans(i,j+1)-vTrans(i,j))       &                   +(vTrans(i,j+1)-vTrans(i,j))*advFac
764       &                   +(rTransKp1(i,j)-rTrans(i,j))*rAdvFac       &                   +(rTransKp1(i,j)-rTrans(i,j))*rAdvFac
765       &                  )*advFac       &                  )*maskInC(i,j,bi,bj)
766       &    )       &    )
767         ENDDO         ENDDO
768        ENDDO        ENDDO
769    
770  #ifdef ALLOW_DEBUG  #ifdef ALLOW_DEBUG
771        IF ( debugLevel .GE. debLevB        IF ( debugLevel .GE. debLevC
772       &   .AND. tracerIdentity.EQ.GAD_TEMPERATURE       &   .AND. trIdentity.EQ.GAD_TEMPERATURE
773       &   .AND. k.EQ.2 .AND. myIter.EQ.1+nIter0       &   .AND. k.EQ.2 .AND. myIter.EQ.1+nIter0
774       &   .AND. nPx.EQ.1 .AND. nPy.EQ.1       &   .AND. nPx.EQ.1 .AND. nPy.EQ.1
775       &   .AND. useCubedSphereExchange ) THEN       &   .AND. useCubedSphereExchange ) THEN
Legend:
Removed from v.1.42  
changed lines
  Added in v.1.61
  ViewVC Help
Powered by ViewVC 1.1.22