/[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.32 by jmc, Thu Dec 16 22:29:34 2004 UTC revision 1.42 by jmc, Sun Nov 19 22:04:09 2006 UTC
# Line 7  CBOP Line 7  CBOP
7  C !ROUTINE: GAD_CALC_RHS  C !ROUTINE: GAD_CALC_RHS
8    
9  C !INTERFACE: ==========================================================  C !INTERFACE: ==========================================================
10        SUBROUTINE GAD_CALC_RHS(        SUBROUTINE GAD_CALC_RHS(
11       I           bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,       I           bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,
12       I           xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,       I           xA, yA, maskUp, uFld, vFld, wFld,
13       I           uVel, vVel, wVel,       I           uTrans, vTrans, rTrans, rTransKp1,
14       I           diffKh, diffK4, KappaR, Tracer,       I           diffKh, diffK4, KappaR, TracerN, TracAB,
15       I           tracerIdentity, advectionScheme, vertAdvecScheme,       I           tracerIdentity, advectionScheme, vertAdvecScheme,
16       I           calcAdvection, implicitAdvection,       I           calcAdvection, implicitAdvection, applyAB_onTracer,
17       U           fVerT, gTracer,       U           fVerT, gTracer,
18       I           myTime, myIter, myThid )       I           myTime, myIter, myThid )
19    
20  C !DESCRIPTION:  C !DESCRIPTION:
21  C Calculates the tendancy of a tracer due to advection and diffusion.  C Calculates the tendency of a tracer due to advection and diffusion.
22  C It calculates the fluxes in each direction indepentently and then  C It calculates the fluxes in each direction indepentently and then
23  C sets the tendancy to the divergence of these fluxes. The advective  C sets the tendency to the divergence of these fluxes. The advective
24  C fluxes are only calculated here when using the linear advection schemes  C fluxes are only calculated here when using the linear advection schemes
25  C otherwise only the diffusive and parameterized fluxes are calculated.  C otherwise only the diffusive and parameterized fluxes are calculated.
26  C  C
# Line 29  C \begin{equation*} Line 29  C \begin{equation*}
29  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}
30  C \end{equation*}  C \end{equation*}
31  C  C
32  C The tendancy is the divergence of the fluxes:  C The tendency is the divergence of the fluxes:
33  C \begin{equation*}  C \begin{equation*}
34  C G_\theta = G_\theta + \nabla \cdot {\bf F}  C G_\theta = G_\theta + \nabla \cdot {\bf F}
35  C \end{equation*}  C \end{equation*}
36  C  C
37  C The tendancy is assumed to contain data on entry.  C The tendency is assumed to contain data on entry.
38    
39  C !USES: ===============================================================  C !USES: ===============================================================
40        IMPLICIT NONE        IMPLICIT NONE
# Line 54  C !INPUT PARAMETERS: =================== Line 54  C !INPUT PARAMETERS: ===================
54  C bi,bj            :: tile indices  C bi,bj            :: tile indices
55  C iMin,iMax        :: loop range for called routines  C iMin,iMax        :: loop range for called routines
56  C jMin,jMax        :: loop range for called routines  C jMin,jMax        :: loop range for called routines
57  C kup              :: index into 2 1/2D array, toggles between 1|2  C k                :: vertical index
58  C kdown            :: index into 2 1/2D array, toggles between 2|1  C kM1              :: =k-1 for k>1, =1 for k=1
59  C kp1              :: =k+1 for k<Nr, =Nr for k=Nr  C kUp              :: index into 2 1/2D array, toggles between 1|2
60    C kDown            :: index into 2 1/2D array, toggles between 2|1
61  C xA,yA            :: areas of X and Y face of tracer cells  C xA,yA            :: areas of X and Y face of tracer cells
62    C maskUp           :: 2-D array for mask at W points
63    C uFld,vFld,wFld   :: Local copy of velocity field (3 components)
64  C uTrans,vTrans    :: 2-D arrays of volume transports at U,V points  C uTrans,vTrans    :: 2-D arrays of volume transports at U,V points
65  C rTrans           :: 2-D arrays of volume transports at W points  C rTrans           :: 2-D arrays of volume transports at W points
66  C rTransKp1        :: 2-D array of volume trans at W pts, interf k+1  C rTransKp1        :: 2-D array of volume trans at W pts, interf k+1
 C maskUp           :: 2-D array for mask at W points  
 C uVel,vVel,wVel   :: 3 components of the velcity field (3-D array)  
67  C diffKh           :: horizontal diffusion coefficient  C diffKh           :: horizontal diffusion coefficient
68  C diffK4           :: bi-harmonic diffusion coefficient  C diffK4           :: bi-harmonic diffusion coefficient
69  C KappaR           :: 2-D array for vertical diffusion coefficient, interf k  C KappaR           :: 2-D array for vertical diffusion coefficient, interf k
70  C Tracer           :: tracer field  C TracerN          :: tracer field @ time-step n (Note: only used
71    C                     if applying AB on tracer field rather than on tendency gTr)
72    C TracAB           :: current tracer field (@ time-step n if applying AB on gTr
73    C                     or extrapolated fwd in time to n+1/2 if applying AB on Tr)
74  C tracerIdentity   :: tracer identifier (required for KPP,GM)  C tracerIdentity   :: tracer identifier (required for KPP,GM)
75  C advectionScheme  :: advection scheme to use (Horizontal plane)  C advectionScheme  :: advection scheme to use (Horizontal plane)
76  C vertAdvecScheme  :: advection scheme to use (Vertical direction)  C vertAdvecScheme  :: advection scheme to use (Vertical direction)
77  C calcAdvection    :: =False if Advec computed with multiDim scheme  C calcAdvection    :: =False if Advec computed with multiDim scheme
78  C implicitAdvection:: =True if vertical Advec computed implicitly  C implicitAdvection:: =True if vertical Advec computed implicitly
79    C applyAB_onTracer :: apply Adams-Bashforth on Tracer (rather than on gTr)
80  C myTime           :: current time  C myTime           :: current time
81  C myIter           :: iteration number  C myIter           :: iteration number
82  C myThid           :: thread number  C myThid           :: thread number
# Line 79  C myThid           :: thread number Line 84  C myThid           :: thread number
84        INTEGER k,kUp,kDown,kM1        INTEGER k,kUp,kDown,kM1
85        _RS xA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS xA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
86        _RS yA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS yA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
87          _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
88          _RL uFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
89          _RL vFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
90          _RL wFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
91        _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
92        _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
93        _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
94        _RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
       _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)  
       _RL uVel  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)  
       _RL vVel  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)  
       _RL wVel  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)  
95        _RL diffKh, diffK4        _RL diffKh, diffK4
96        _RL KappaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL KappaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
97        _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)
98          _RL TracAB (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
99        INTEGER tracerIdentity        INTEGER tracerIdentity
100        INTEGER advectionScheme, vertAdvecScheme        INTEGER advectionScheme, vertAdvecScheme
101        LOGICAL calcAdvection        LOGICAL calcAdvection
102        LOGICAL implicitAdvection        LOGICAL implicitAdvection, applyAB_onTracer
103        _RL     myTime        _RL     myTime
104        INTEGER myIter, myThid        INTEGER myIter, myThid
105    
106  C !OUTPUT PARAMETERS: ==================================================  C !OUTPUT PARAMETERS: ==================================================
107  C gTracer          :: tendancy array  C gTracer          :: tendency array
108  C fVerT            :: 2 1/2D arrays for vertical advective flux  C fVerT            :: 2 1/2D arrays for vertical advective flux
109        _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)
110        _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
# Line 111  C fMer             :: meridional flux Line 117  C fMer             :: meridional flux
117  C af               :: advective flux  C af               :: advective flux
118  C df               :: diffusive flux  C df               :: diffusive flux
119  C localT           :: local copy of tracer field  C localT           :: local copy of tracer field
120    C locABT           :: local copy of (AB-extrapolated) tracer field
121  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
       INTEGER     kk  
122        CHARACTER*8 diagName        CHARACTER*8 diagName
123        CHARACTER*4 GAD_DIAG_SUFX, diagSufx        CHARACTER*4 GAD_DIAG_SUFX, diagSufx
124        EXTERNAL    GAD_DIAG_SUFX        EXTERNAL    GAD_DIAG_SUFX
125  #endif  #endif
126        INTEGER i,j        INTEGER i,j
# Line 124  C localT           :: local copy of trac Line 130  C localT           :: local copy of trac
130        _RL af    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL af    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
131        _RL df    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL df    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
132        _RL localT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL localT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
133          _RL locABT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
134        _RL advFac, rAdvFac        _RL advFac, rAdvFac
135  CEOP  CEOP
136    
# Line 142  C--   Set diagnostic suffix for the curr Line 149  C--   Set diagnostic suffix for the curr
149    
150        advFac  = 0. _d 0        advFac  = 0. _d 0
151        IF (calcAdvection) advFac = 1. _d 0        IF (calcAdvection) advFac = 1. _d 0
152        rAdvFac = rkFac*advFac        rAdvFac = rkSign*advFac
153        IF (implicitAdvection) rAdvFac = 0. _d 0        IF (implicitAdvection) rAdvFac = 0. _d 0
154    
155        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
# Line 156  C--   Set diagnostic suffix for the curr Line 163  C--   Set diagnostic suffix for the curr
163        ENDDO        ENDDO
164    
165  C--   Make local copy of tracer array  C--   Make local copy of tracer array
166        DO j=1-OLy,sNy+OLy        IF ( applyAB_onTracer ) THEN
167         DO i=1-OLx,sNx+OLx          DO j=1-OLy,sNy+OLy
168          localT(i,j)=tracer(i,j,k,bi,bj)           DO i=1-OLx,sNx+OLx
169         ENDDO            localT(i,j)=TracerN(i,j,k,bi,bj)
170        ENDDO            locABT(i,j)= TracAB(i,j,k,bi,bj)
171             ENDDO
172            ENDDO
173          ELSE
174            DO j=1-OLy,sNy+OLy
175             DO i=1-OLx,sNx+OLx
176              localT(i,j)= TracAB(i,j,k,bi,bj)
177              locABT(i,j)= TracAB(i,j,k,bi,bj)
178             ENDDO
179            ENDDO
180          ENDIF
181    
182  C--   Unless we have already calculated the advection terms we initialize  C--   Unless we have already calculated the advection terms we initialize
183  C     the tendency to zero.  C     the tendency to zero.
# Line 190  C--   Initialize net flux in X direction Line 207  C--   Initialize net flux in X direction
207  C-    Advective flux in X  C-    Advective flux in X
208        IF (calcAdvection) THEN        IF (calcAdvection) THEN
209          IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN          IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
210            CALL GAD_C2_ADV_X(bi,bj,k,uTrans,localT,af,myThid)            CALL GAD_C2_ADV_X(bi,bj,k,uTrans,locABT,af,myThid)
211            ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST
212         &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
213              CALL GAD_DST2U1_ADV_X( bi,bj,k, advectionScheme,
214         I            dTtracerLev(k), uTrans, uFld, locABT,
215         O            af, myThid )
216          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
217            CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, dTtracerLev(k),            CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, dTtracerLev(k),
218       I            uTrans, uVel, maskW(1-Olx,1-Oly,k,bi,bj), localT,       I            uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,
219       O            af, myThid )       O            af, myThid )
220          ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN          ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
221            CALL GAD_U3_ADV_X(bi,bj,k,uTrans,localT,af,myThid)            CALL GAD_U3_ADV_X(bi,bj,k,uTrans,locABT,af,myThid)
222          ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN          ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
223            CALL GAD_C4_ADV_X(bi,bj,k,uTrans,localT,af,myThid)            CALL GAD_C4_ADV_X(bi,bj,k,uTrans,locABT,af,myThid)
224          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
225            CALL GAD_DST3_ADV_X( bi,bj,k, dTtracerLev(k),            CALL GAD_DST3_ADV_X( bi,bj,k, dTtracerLev(k),
226       I            uTrans, uVel, maskW(1-Olx,1-Oly,k,bi,bj), localT,       I            uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,
227       O            af, myThid )       O            af, myThid )
228          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
229             IF ( inAdMode ) THEN
230    cph This block is to trick the adjoint:
231    cph IF inAdExact=.FALSE., we want to use DST3
232    cph with limiters in forward, but without limiters in reverse.
233              CALL GAD_DST3_ADV_X( bi,bj,k, dTtracerLev(k),
234         I           uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,
235         O           af, myThid )
236             ELSE
237            CALL GAD_DST3FL_ADV_X( bi,bj,k, dTtracerLev(k),            CALL GAD_DST3FL_ADV_X( bi,bj,k, dTtracerLev(k),
238       I            uTrans, uVel, maskW(1-Olx,1-Oly,k,bi,bj), localT,       I           uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,
239       O            af, myThid )       O           af, myThid )
240             ENDIF
241          ELSE          ELSE
242           STOP 'GAD_CALC_RHS: Bad advectionScheme (X)'           STOP 'GAD_CALC_RHS: Bad advectionScheme (X)'
243          ENDIF          ENDIF
# Line 218  C-    Advective flux in X Line 249  C-    Advective flux in X
249  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
250          IF ( useDiagnostics ) THEN          IF ( useDiagnostics ) THEN
251            diagName = 'ADVx'//diagSufx            diagName = 'ADVx'//diagSufx
252            kk = -k            CALL DIAGNOSTICS_FILL(af,diagName, k,1, 2,bi,bj, myThid)
           CALL DIAGNOSTICS_FILL(af,diagName, kk,1, 2,bi,bj, myThid)  
253          ENDIF          ENDIF
254  #endif  #endif
255        ENDIF        ENDIF
# Line 243  C-    Add bi-harmonic diffusive flux in Line 273  C-    Add bi-harmonic diffusive flux in
273  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
274  C-    GM/Redi flux in X  C-    GM/Redi flux in X
275        IF (useGMRedi) THEN        IF (useGMRedi) THEN
276  C *note* should update GMREDI_XTRANSPORT to use localT and set df  *aja*  C *note* should update GMREDI_XTRANSPORT to set df  *aja*
277          CALL GMREDI_XTRANSPORT(          IF ( applyAB_onTracer ) THEN
278       I     iMin,iMax,jMin,jMax,bi,bj,K,            CALL GMREDI_XTRANSPORT(
279       I     xA,Tracer,tracerIdentity,       I         iMin,iMax,jMin,jMax,bi,bj,k,
280       U     df,       I         xA,TracerN,tracerIdentity,
281       I     myThid)       U         df,
282         I         myThid)
283            ELSE
284              CALL GMREDI_XTRANSPORT(
285         I         iMin,iMax,jMin,jMax,bi,bj,k,
286         I         xA,TracAB, tracerIdentity,
287         U         df,
288         I         myThid)
289            ENDIF
290        ENDIF        ENDIF
291  #endif  #endif
292        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
# Line 262  C-    Diagnostics of Tracer flux in X di Line 300  C-    Diagnostics of Tracer flux in X di
300  C       excluding advective terms:  C       excluding advective terms:
301        IF ( useDiagnostics .AND.        IF ( useDiagnostics .AND.
302       &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. useGMRedi) ) THEN       &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. useGMRedi) ) THEN
303            diagName = 'DIFx'//diagSufx            diagName = 'DFxE'//diagSufx
304            kk = -k            CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
           CALL DIAGNOSTICS_FILL(df,diagName, kk,1, 2,bi,bj, myThid)  
305        ENDIF        ENDIF
306  #endif  #endif
307    
# Line 278  C--   Initialize net flux in Y direction Line 315  C--   Initialize net flux in Y direction
315  C-    Advective flux in Y  C-    Advective flux in Y
316        IF (calcAdvection) THEN        IF (calcAdvection) THEN
317          IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN          IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
318            CALL GAD_C2_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)            CALL GAD_C2_ADV_Y(bi,bj,k,vTrans,locABT,af,myThid)
319            ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST
320         &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
321              CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme,
322         I            dTtracerLev(k), vTrans, vFld, locABT,
323         O            af, myThid )
324          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
325            CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, dTtracerLev(k),            CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, dTtracerLev(k),
326       I            vTrans, vVel, maskS(1-Olx,1-Oly,k,bi,bj), localT,       I            vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,
327       O            af, myThid )       O            af, myThid )
328          ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN          ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
329            CALL GAD_U3_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)            CALL GAD_U3_ADV_Y(bi,bj,k,vTrans,locABT,af,myThid)
330          ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN          ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
331            CALL GAD_C4_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)            CALL GAD_C4_ADV_Y(bi,bj,k,vTrans,locABT,af,myThid)
332          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
333            CALL GAD_DST3_ADV_Y( bi,bj,k, dTtracerLev(k),            CALL GAD_DST3_ADV_Y( bi,bj,k, dTtracerLev(k),
334       I            vTrans, vVel, maskS(1-Olx,1-Oly,k,bi,bj), localT,       I            vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,
335       O            af, myThid )       O            af, myThid )
336          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
337             IF ( inAdMode ) THEN
338    cph This block is to trick the adjoint:
339    cph IF inAdExact=.FALSE., we want to use DST3
340    cph with limiters in forward, but without limiters in reverse.
341              CALL GAD_DST3_ADV_Y( bi,bj,k, dTtracerLev(k),
342         I           vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,
343         O           af, myThid )
344             ELSE
345            CALL GAD_DST3FL_ADV_Y( bi,bj,k, dTtracerLev(k),            CALL GAD_DST3FL_ADV_Y( bi,bj,k, dTtracerLev(k),
346       I            vTrans, vVel, maskS(1-Olx,1-Oly,k,bi,bj), localT,       I           vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,
347       O            af, myThid )       O           af, myThid )
348             ENDIF
349          ELSE          ELSE
350            STOP 'GAD_CALC_RHS: Bad advectionScheme (Y)'            STOP 'GAD_CALC_RHS: Bad advectionScheme (Y)'
351          ENDIF          ENDIF
# Line 306  C-    Advective flux in Y Line 357  C-    Advective flux in Y
357  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
358          IF ( useDiagnostics ) THEN          IF ( useDiagnostics ) THEN
359            diagName = 'ADVy'//diagSufx            diagName = 'ADVy'//diagSufx
360            kk = -k            CALL DIAGNOSTICS_FILL(af,diagName, k,1, 2,bi,bj, myThid)
           CALL DIAGNOSTICS_FILL(af,diagName, kk,1, 2,bi,bj, myThid)  
361          ENDIF          ENDIF
362  #endif  #endif
363        ENDIF        ENDIF
# Line 331  C-    Add bi-harmonic flux in Y Line 381  C-    Add bi-harmonic flux in Y
381  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
382  C-    GM/Redi flux in Y  C-    GM/Redi flux in Y
383        IF (useGMRedi) THEN        IF (useGMRedi) THEN
384  C *note* should update GMREDI_YTRANSPORT to use localT and set df  *aja*  C *note* should update GMREDI_YTRANSPORT to set df  *aja*
385         CALL GMREDI_YTRANSPORT(          IF ( applyAB_onTracer ) THEN
386       I     iMin,iMax,jMin,jMax,bi,bj,K,            CALL GMREDI_YTRANSPORT(
387       I     yA,Tracer,tracerIdentity,       I         iMin,iMax,jMin,jMax,bi,bj,k,
388       U     df,       I         yA,TracerN,tracerIdentity,
389       I     myThid)       U         df,
390         I         myThid)
391            ELSE
392              CALL GMREDI_YTRANSPORT(
393         I         iMin,iMax,jMin,jMax,bi,bj,k,
394         I         yA,TracAB, tracerIdentity,
395         U         df,
396         I         myThid)
397            ENDIF
398        ENDIF        ENDIF
399  #endif  #endif
400        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
# Line 350  C-    Diagnostics of Tracer flux in Y di Line 408  C-    Diagnostics of Tracer flux in Y di
408  C       excluding advective terms:  C       excluding advective terms:
409        IF ( useDiagnostics .AND.        IF ( useDiagnostics .AND.
410       &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. useGMRedi) ) THEN       &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. useGMRedi) ) THEN
411            diagName = 'DIFy'//diagSufx            diagName = 'DFyE'//diagSufx
412            kk = -k            CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
           CALL DIAGNOSTICS_FILL(df,diagName, kk,1, 2,bi,bj, myThid)  
413        ENDIF        ENDIF
414  #endif  #endif
415    
# Line 360  C--   Compute vertical flux fVerT(kUp) a Line 417  C--   Compute vertical flux fVerT(kUp) a
417  C-    Advective flux in R  C-    Advective flux in R
418  #ifdef ALLOW_AIM  #ifdef ALLOW_AIM
419  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
420        IF (calcAdvection .AND. .NOT.implicitAdvection .AND. K.GE.2 .AND.        IF (calcAdvection .AND. .NOT.implicitAdvection .AND. k.GE.2 .AND.
421       &     (.NOT.useAIM .OR.tracerIdentity.NE.GAD_SALINITY .OR.K.LT.Nr)       &     (.NOT.useAIM .OR.tracerIdentity.NE.GAD_SALINITY .OR.k.LT.Nr)
422       &   ) THEN       &   ) THEN
423  #else  #else
424        IF (calcAdvection .AND. .NOT.implicitAdvection .AND. K.GE.2) THEN        IF (calcAdvection .AND. .NOT.implicitAdvection .AND. k.GE.2) THEN
425  #endif  #endif
426  C-    Compute vertical advective flux in the interior:  C-    Compute vertical advective flux in the interior:
427          IF (vertAdvecScheme.EQ.ENUM_CENTERED_2ND) THEN          IF (vertAdvecScheme.EQ.ENUM_CENTERED_2ND) THEN
428            CALL GAD_C2_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)            CALL GAD_C2_ADV_R(bi,bj,k,rTrans,TracAB,af,myThid)
429            ELSEIF ( vertAdvecScheme.EQ.ENUM_UPWIND_1RST
430         &          .OR. vertAdvecScheme.EQ.ENUM_DST2 ) THEN
431              CALL GAD_DST2U1_ADV_R( bi,bj,k, vertAdvecScheme,
432         I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),
433         O         af, myThid )
434          ELSEIF (vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF (vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN
435            CALL GAD_FLUXLIMIT_ADV_R(            CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k,
436       &         bi,bj,k,dTtracerLev(k),rTrans,wVel,tracer,af,myThid)       I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),
437         O         af, myThid )
438          ELSEIF (vertAdvecScheme.EQ.ENUM_UPWIND_3RD ) THEN          ELSEIF (vertAdvecScheme.EQ.ENUM_UPWIND_3RD ) THEN
439            CALL GAD_U3_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)            CALL GAD_U3_ADV_R(bi,bj,k,rTrans,TracAB,af,myThid)
440          ELSEIF (vertAdvecScheme.EQ.ENUM_CENTERED_4TH) THEN          ELSEIF (vertAdvecScheme.EQ.ENUM_CENTERED_4TH) THEN
441            CALL GAD_C4_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)            CALL GAD_C4_ADV_R(bi,bj,k,rTrans,TracAB,af,myThid)
442          ELSEIF (vertAdvecScheme.EQ.ENUM_DST3 ) THEN          ELSEIF (vertAdvecScheme.EQ.ENUM_DST3 ) THEN
443            CALL GAD_DST3_ADV_R(            CALL GAD_DST3_ADV_R( bi,bj,k,
444       &         bi,bj,k,dTtracerLev(k),rTrans,wVel,tracer,af,myThid)       I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),
445         O         af, myThid )
446          ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN          ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
447            CALL GAD_DST3FL_ADV_R(  cph This block is to trick the adjoint:
448       &         bi,bj,k,dTtracerLev(k),rTrans,wVel,tracer,af,myThid)  cph IF inAdExact=.FALSE., we want to use DST3
449    cph with limiters in forward, but without limiters in reverse.
450              IF ( inAdMode ) THEN
451               CALL GAD_DST3_ADV_R( bi,bj,k,
452         I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),
453         O         af, myThid )
454              ELSE
455               CALL GAD_DST3FL_ADV_R( bi,bj,k,
456         I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),
457         O         af, myThid )
458              ENDIF
459          ELSE          ELSE
460            STOP 'GAD_CALC_RHS: Bad vertAdvecScheme (R)'            STOP 'GAD_CALC_RHS: Bad vertAdvecScheme (R)'
461          ENDIF          ENDIF
# Line 394  C-     add the advective flux to fVerT Line 468  C-     add the advective flux to fVerT
468  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
469          IF ( useDiagnostics ) THEN          IF ( useDiagnostics ) THEN
470            diagName = 'ADVr'//diagSufx            diagName = 'ADVr'//diagSufx
471            kk = -k            CALL DIAGNOSTICS_FILL(af,diagName, k,1, 2,bi,bj, myThid)
472            CALL DIAGNOSTICS_FILL(af,diagName, kk,1, 2,bi,bj, myThid)  C- note: needs to explicitly increment the counter since DIAGNOSTICS_FILL
473    C        does it only if k=1 (never the case here)
474              IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT(diagName,bi,bj,myThid)
475          ENDIF          ENDIF
476  #endif  #endif
477        ENDIF        ENDIF
# Line 410  C           boundary condition. Line 486  C           boundary condition.
486          ENDDO          ENDDO
487         ENDDO         ENDDO
488        ELSE        ELSE
489         CALL GAD_DIFF_R(bi,bj,k,KappaR,tracer,df,myThid)         IF ( applyAB_onTracer ) THEN
490             CALL GAD_DIFF_R(bi,bj,k,KappaR,TracerN,df,myThid)
491           ELSE
492             CALL GAD_DIFF_R(bi,bj,k,KappaR,TracAB, df,myThid)
493           ENDIF
494        ENDIF        ENDIF
495    
496  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
497  C-    GM/Redi flux in R  C-    GM/Redi flux in R
498        IF (useGMRedi) THEN        IF (useGMRedi) THEN
499  C *note* should update GMREDI_RTRANSPORT to set df  *aja*  C *note* should update GMREDI_RTRANSPORT to set df  *aja*
500         CALL GMREDI_RTRANSPORT(          IF ( applyAB_onTracer ) THEN
501       I     iMin,iMax,jMin,jMax,bi,bj,K,            CALL GMREDI_RTRANSPORT(
502       I     Tracer,tracerIdentity,       I         iMin,iMax,jMin,jMax,bi,bj,k,
503       U     df,       I         TracerN,tracerIdentity,
504       I     myThid)       U         df,
505         I         myThid)
506            ELSE
507              CALL GMREDI_RTRANSPORT(
508         I         iMin,iMax,jMin,jMax,bi,bj,k,
509         I         TracAB, tracerIdentity,
510         U         df,
511         I         myThid)
512            ENDIF
513        ENDIF        ENDIF
514  #endif  #endif
515    
# Line 432  C *note* should update GMREDI_RTRANSPORT Line 520  C *note* should update GMREDI_RTRANSPORT
520        ENDDO        ENDDO
521    
522  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
523  C-    Diagnostics of Tracer flux in R dir (mainly Diffusive terms),  C-    Diagnostics of Tracer flux in R dir (mainly Diffusive terms),
524  C       Explicit terms only & excluding advective terms:  C       Explicit terms only & excluding advective terms:
525        IF ( useDiagnostics .AND.        IF ( useDiagnostics .AND.
526       &    (.NOT.implicitDiffusion .OR. useGMRedi) ) THEN       &    (.NOT.implicitDiffusion .OR. useGMRedi) ) THEN
527            diagName = 'DFrE'//diagSufx            diagName = 'DFrE'//diagSufx
528            kk = -k            CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
           CALL DIAGNOSTICS_FILL(df,diagName, kk,1, 2,bi,bj, myThid)  
529        ENDIF        ENDIF
530  #endif  #endif
531    
# Line 484  C--   Divergence of fluxes Line 571  C--   Divergence of fluxes
571       &   -_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)*recip_rA(i,j,bi,bj)
572       &   *( (fZon(i+1,j)-fZon(i,j))       &   *( (fZon(i+1,j)-fZon(i,j))
573       &     +(fMer(i,j+1)-fMer(i,j))       &     +(fMer(i,j+1)-fMer(i,j))
574       &     +(fVerT(i,j,kUp)-fVerT(i,j,kDown))*rkFac       &     +(fVerT(i,j,kDown)-fVerT(i,j,kUp))*rkSign
575       &     -localT(i,j)*( (uTrans(i+1,j)-uTrans(i,j))       &     -localT(i,j)*( (uTrans(i+1,j)-uTrans(i,j))
576       &                   +(vTrans(i,j+1)-vTrans(i,j))       &                   +(vTrans(i,j+1)-vTrans(i,j))
577       &                   +(rTrans(i,j)-rTransKp1(i,j))*rAdvFac       &                   +(rTransKp1(i,j)-rTrans(i,j))*rAdvFac
578       &                  )*advFac       &                  )*advFac
579       &    )       &    )
580         ENDDO         ENDDO
# Line 503  C--   Divergence of fluxes Line 590  C--   Divergence of fluxes
590       &             fZon,fMer, k, standardMessageUnit,bi,bj,myThid )       &             fZon,fMer, k, standardMessageUnit,bi,bj,myThid )
591        ENDIF        ENDIF
592  #endif /* ALLOW_DEBUG */  #endif /* ALLOW_DEBUG */
593    
594        RETURN        RETURN
595        END        END
Legend:
Removed from v.1.32  
changed lines
  Added in v.1.42
  ViewVC Help
Powered by ViewVC 1.1.22