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