/[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.6 by adcroft, Tue Sep 4 17:16:11 2001 UTC revision 1.44 by adcroft, Sat Jan 20 21:20:11 2007 UTC
# Line 3  C $Name$ Line 3  C $Name$
3    
4  #include "GAD_OPTIONS.h"  #include "GAD_OPTIONS.h"
5    
6        SUBROUTINE GAD_CALC_RHS(  CBOP
7    C !ROUTINE: GAD_CALC_RHS
8    
9    C !INTERFACE: ==========================================================
10          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,       I           diffKh, diffK4, KappaR, TracerN, TracAB,
15         I           tracerIdentity, advectionScheme, vertAdvecScheme,
16         I           calcAdvection, implicitAdvection, applyAB_onTracer,
17       U           fVerT, gTracer,       U           fVerT, gTracer,
18       I           myThid )       I           myTime, myIter, myThid )
 C     /==========================================================\  
 C     | SUBROUTINE GAD_CALC_RHS                                  |  
 C     |==========================================================|  
 C     \==========================================================/  
       IMPLICIT NONE  
19    
20  C     == GLobal variables ==  C !DESCRIPTION:
21    C Calculates the tendency of a tracer due to advection and diffusion.
22    C It calculates the fluxes in each direction indepentently and then
23    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
25    C otherwise only the diffusive and parameterized fluxes are calculated.
26    C
27    C Contributions to the flux are calculated and added:
28    C \begin{equation*}
29    C {\bf F} = {\bf F}_{adv} + {\bf F}_{diff} +{\bf F}_{GM} + {\bf F}_{KPP}
30    C \end{equation*}
31    C
32    C The tendency is the divergence of the fluxes:
33    C \begin{equation*}
34    C G_\theta = G_\theta + \nabla \cdot {\bf F}
35    C \end{equation*}
36    C
37    C The tendency is assumed to contain data on entry.
38    
39    C !USES: ===============================================================
40          IMPLICIT NONE
41  #include "SIZE.h"  #include "SIZE.h"
42  #include "EEPARAMS.h"  #include "EEPARAMS.h"
43  #include "PARAMS.h"  #include "PARAMS.h"
44  #include "GRID.h"  #include "GRID.h"
45  #include "DYNVARS.h"  #include "SURFACE.h"
46  #include "GAD.h"  #include "GAD.h"
47    
48  C     == Routine arguments ==  #ifdef ALLOW_AUTODIFF_TAMC
49        INTEGER k,kUp,kDown,kM1  #include "tamc.h"
50    #include "tamc_keys.h"
51    #endif /* ALLOW_AUTODIFF_TAMC */
52    
53    C !INPUT PARAMETERS: ===================================================
54    C bi,bj            :: tile indices
55    C iMin,iMax        :: loop range for called routines
56    C jMin,jMax        :: loop range for called routines
57    C k                :: vertical index
58    C kM1              :: =k-1 for k>1, =1 for k=1
59    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
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
65    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
67    C diffKh           :: horizontal diffusion coefficient
68    C diffK4           :: bi-harmonic diffusion coefficient
69    C KappaR           :: 2-D array for vertical diffusion coefficient, interf k
70    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)
75    C advectionScheme  :: advection scheme to use (Horizontal plane)
76    C vertAdvecScheme  :: advection scheme to use (Vertical direction)
77    C calcAdvection    :: =False if Advec computed with multiDim scheme
78    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
81    C myIter           :: iteration number
82    C myThid           :: thread number
83        INTEGER bi,bj,iMin,iMax,jMin,jMax        INTEGER bi,bj,iMin,iMax,jMin,jMax
84          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        _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
95        _RL diffKh, diffK4        _RL diffKh, diffK4
96        _RL KappaRT(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _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        INTEGER advectionScheme, vertAdvecScheme
101        _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        LOGICAL calcAdvection
102          LOGICAL implicitAdvection, applyAB_onTracer
103          _RL     myTime
104          INTEGER myIter, myThid
105    
106    C !OUTPUT PARAMETERS: ==================================================
107    C gTracer          :: tendency array
108    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        INTEGER myThid        _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
111    
112  C     == Local variables ==  C !LOCAL VARIABLES: ====================================================
113  C     I, J, K - Loop counters  C i,j              :: loop indices
114    C df4              :: used for storing del^2 T for bi-harmonic term
115    C fZon             :: zonal flux
116    C fMer             :: meridional flux
117    C af               :: advective flux
118    C df               :: diffusive flux
119    C localT           :: local copy of tracer field
120    C locABT           :: local copy of (AB-extrapolated) tracer field
121    #ifdef ALLOW_DIAGNOSTICS
122          CHARACTER*8 diagName
123          CHARACTER*4 GAD_DIAG_SUFX, diagSufx
124          EXTERNAL    GAD_DIAG_SUFX
125    #endif
126        INTEGER i,j        INTEGER i,j
       LOGICAL TOP_LAYER  
       _RL afFacT, dfFacT  
127        _RL df4   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL df4   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
128        _RL fZon  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fZon  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
129        _RL fMer  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fMer  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
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
135    CEOP
136    
137  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
138  C--   only the kUp part of fverT is set in this subroutine  C--   only the kUp part of fverT is set in this subroutine
139  C--   the kDown is still required  C--   the kDown is still required
140        fVerT(1,1,kDown) = fVerT(1,1,kDown)        fVerT(1,1,kDown) = fVerT(1,1,kDown)
141  #endif  #endif
       DO j=1-OLy,sNy+OLy  
        DO i=1-OLx,sNx+OLx  
         fZon(i,j)      = 0.0  
         fMer(i,j)      = 0.0  
         fVerT(i,j,kUp) = 0.0  
        ENDDO  
       ENDDO  
142    
143        afFacT = 1. _d 0  #ifdef ALLOW_DIAGNOSTICS
144        dfFacT = 1. _d 0  C--   Set diagnostic suffix for the current tracer
145        TOP_LAYER = K .EQ. 1        IF ( useDiagnostics ) THEN
146            diagSufx = GAD_DIAG_SUFX( tracerIdentity, myThid )
147          ENDIF
148    #endif
149    
150          advFac  = 0. _d 0
151          IF (calcAdvection) advFac = 1. _d 0
152          rAdvFac = rkSign*advFac
153          IF (implicitAdvection) rAdvFac = 0. _d 0
154    
 C--   Make local copy of tracer array  
155        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
156         DO i=1-OLx,sNx+OLx         DO i=1-OLx,sNx+OLx
157          localT(i,j)=tracer(i,j,k,bi,bj)          fZon(i,j)      = 0. _d 0
158            fMer(i,j)      = 0. _d 0
159            fVerT(i,j,kUp) = 0. _d 0
160            df(i,j)        = 0. _d 0
161            df4(i,j)       = 0. _d 0
162         ENDDO         ENDDO
163        ENDDO        ENDDO
164    
165    C--   Make local copy of tracer array
166          IF ( applyAB_onTracer ) THEN
167            DO j=1-OLy,sNy+OLy
168             DO i=1-OLx,sNx+OLx
169              localT(i,j)=TracerN(i,j,k,bi,bj)
170              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
183    C     the tendency to zero.
184    C     <== now done earlier at the beginning of thermodynamics.
185    c     IF (calcAdvection) THEN
186    c      DO j=1-Oly,sNy+Oly
187    c       DO i=1-Olx,sNx+Olx
188    c        gTracer(i,j,k,bi,bj)=0. _d 0
189    c       ENDDO
190    c      ENDDO
191    c     ENDIF
192    
193  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
194        IF (diffK4 .NE. 0.) THEN        IF (diffK4 .NE. 0.) THEN
# Line 89  C--   Pre-calculate del^2 T if bi-harmon Line 200  C--   Pre-calculate del^2 T if bi-harmon
200  C--   Initialize net flux in X direction  C--   Initialize net flux in X direction
201        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
202         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
203          fZon(i,j) = 0.          fZon(i,j) = 0. _d 0
204         ENDDO         ENDDO
205        ENDDO        ENDDO
206    
207  C-    Advective flux in X  C-    Advective flux in X
208        IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN        IF (calcAdvection) THEN
209         CALL GAD_C2_ADV_X(bi,bj,k,uTrans,localT,af,myThid)          IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
210        ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN            CALL GAD_C2_ADV_X(bi,bj,k,uTrans,locABT,af,myThid)
211         CALL GAD_FLUXLIMIT_ADV_X(          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST
212       &      bi,bj,k,deltaTtracer,uTrans,uVel,localT,af,myThid)       &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
213        ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN            CALL GAD_DST2U1_ADV_X( bi,bj,k, advectionScheme,
214         CALL GAD_U3_ADV_X(bi,bj,k,uTrans,localT,af,myThid)       I            dTtracerLev(k), uTrans, uFld, locABT,
215        ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN       O            af, myThid )
216         CALL GAD_C4_ADV_X(bi,bj,k,uTrans,localT,af,myThid)          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
217        ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN            CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, dTtracerLev(k),
218         CALL GAD_DST3_ADV_X(       I            uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,
219       &       bi,bj,k,deltaTtracer,uTrans,uVel,localT,af,myThid)       O            af, myThid )
220        ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN          ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
221         CALL GAD_DST3FL_ADV_X(            CALL GAD_U3_ADV_X(bi,bj,k,uTrans,locABT,af,myThid)
222       &       bi,bj,k,deltaTtracer,uTrans,uVel,localT,af,myThid)          ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
223        ELSE            CALL GAD_C4_ADV_X(bi,bj,k,uTrans,locABT,af,myThid)
224         STOP 'GAD_CALC_RHS: Bad advectionScheme (X)'          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
225              CALL GAD_DST3_ADV_X( bi,bj,k, dTtracerLev(k),
226         I            uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,
227         O            af, myThid )
228            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),
238         I           uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,
239         O           af, myThid )
240             ENDIF
241            ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN
242              CALL GAD_OS7MP_ADV_X( bi,bj,k, dTtracerLev(k),
243         I            uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,
244         O            af, myThid )
245            ELSE
246             STOP 'GAD_CALC_RHS: Bad advectionScheme (X)'
247            ENDIF
248            DO j=1-Oly,sNy+Oly
249             DO i=1-Olx,sNx+Olx
250              fZon(i,j) = fZon(i,j) + af(i,j)
251             ENDDO
252            ENDDO
253    #ifdef ALLOW_DIAGNOSTICS
254            IF ( useDiagnostics ) THEN
255              diagName = 'ADVx'//diagSufx
256              CALL DIAGNOSTICS_FILL(af,diagName, k,1, 2,bi,bj, myThid)
257            ENDIF
258    #endif
259        ENDIF        ENDIF
       DO j=1-Oly,sNy+Oly  
        DO i=1-Olx,sNx+Olx  
         fZon(i,j) = fZon(i,j) + af(i,j)  
        ENDDO  
       ENDDO  
260    
261  C-    Diffusive flux in X  C-    Diffusive flux in X
262        IF (diffKh.NE.0.) THEN        IF (diffKh.NE.0.) THEN
# Line 124  C-    Diffusive flux in X Line 264  C-    Diffusive flux in X
264        ELSE        ELSE
265         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
266          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
267           df(i,j) = 0.           df(i,j) = 0. _d 0
268          ENDDO          ENDDO
269         ENDDO         ENDDO
270        ENDIF        ENDIF
271    
272    C-    Add bi-harmonic diffusive flux in X
273          IF (diffK4 .NE. 0.) THEN
274           CALL GAD_BIHARM_X(bi,bj,k,xA,df4,diffK4,df,myThid)
275          ENDIF
276    
277  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
278  C-    GM/Redi flux in X  C-    GM/Redi flux in X
279        IF (useGMRedi) THEN        IF (useGMRedi) THEN
280  C *note* should update GMREDI_XTRANSPORT to use localT and set df  *aja*  C *note* should update GMREDI_XTRANSPORT to set df  *aja*
281          CALL GMREDI_XTRANSPORT(          IF ( applyAB_onTracer ) THEN
282       I     iMin,iMax,jMin,jMax,bi,bj,K,            CALL GMREDI_XTRANSPORT(
283       I     xA,Tracer,       I         iMin,iMax,jMin,jMax,bi,bj,k,
284       U     df,       I         xA,TracerN,tracerIdentity,
285       I     myThid)       U         df,
286         I         myThid)
287            ELSE
288              CALL GMREDI_XTRANSPORT(
289         I         iMin,iMax,jMin,jMax,bi,bj,k,
290         I         xA,TracAB, tracerIdentity,
291         U         df,
292         I         myThid)
293            ENDIF
294        ENDIF        ENDIF
295  #endif  #endif
296    C     anelastic: advect.fluxes are scaled by rhoFac but hor.diff. flx are not
297        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
298         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
299          fZon(i,j) = fZon(i,j) + df(i,j)          fZon(i,j) = fZon(i,j) + df(i,j)*rhoFacC(k)
300         ENDDO         ENDDO
301        ENDDO        ENDDO
302    
303  C-    Bi-harmonic duffusive flux in X  #ifdef ALLOW_DIAGNOSTICS
304        IF (diffK4 .NE. 0.) THEN  C-    Diagnostics of Tracer flux in X dir (mainly Diffusive term),
305         CALL GAD_BIHARM_X(bi,bj,k,xA,df4,diffK4,df,myThid)  C       excluding advective terms:
306         DO j=1-Oly,sNy+Oly        IF ( useDiagnostics .AND.
307          DO i=1-Olx,sNx+Olx       &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. useGMRedi) ) THEN
308           fZon(i,j) = fZon(i,j) + df(i,j)            diagName = 'DFxE'//diagSufx
309          ENDDO            CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
        ENDDO  
310        ENDIF        ENDIF
311    #endif
312    
313  C--   Initialize net flux in Y direction  C--   Initialize net flux in Y direction
314        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
315         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
316          fMer(i,j) = 0.          fMer(i,j) = 0. _d 0
317         ENDDO         ENDDO
318        ENDDO        ENDDO
319    
320  C-    Advective flux in Y  C-    Advective flux in Y
321        IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN        IF (calcAdvection) THEN
322         CALL GAD_C2_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)          IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
323        ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN            CALL GAD_C2_ADV_Y(bi,bj,k,vTrans,locABT,af,myThid)
324         CALL GAD_FLUXLIMIT_ADV_Y(          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST
325       &       bi,bj,k,deltaTtracer,vTrans,vVel,localT,af,myThid)       &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
326        ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN            CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme,
327         CALL GAD_U3_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)       I            dTtracerLev(k), vTrans, vFld, locABT,
328        ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN       O            af, myThid )
329         CALL GAD_C4_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
330        ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN            CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, dTtracerLev(k),
331         CALL GAD_DST3_ADV_Y(       I            vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,
332       &       bi,bj,k,deltaTtracer,vTrans,vVel,localT,af,myThid)       O            af, myThid )
333        ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN          ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
334         CALL GAD_DST3FL_ADV_Y(            CALL GAD_U3_ADV_Y(bi,bj,k,vTrans,locABT,af,myThid)
335       &       bi,bj,k,deltaTtracer,vTrans,vVel,localT,af,myThid)          ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
336        ELSE            CALL GAD_C4_ADV_Y(bi,bj,k,vTrans,locABT,af,myThid)
337         STOP 'GAD_CALC_RHS: Bad advectionScheme (Y)'          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
338              CALL GAD_DST3_ADV_Y( bi,bj,k, dTtracerLev(k),
339         I            vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,
340         O            af, myThid )
341            ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
342             IF ( inAdMode ) THEN
343    cph This block is to trick the adjoint:
344    cph IF inAdExact=.FALSE., we want to use DST3
345    cph with limiters in forward, but without limiters in reverse.
346              CALL GAD_DST3_ADV_Y( bi,bj,k, dTtracerLev(k),
347         I           vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,
348         O           af, myThid )
349             ELSE
350              CALL GAD_DST3FL_ADV_Y( bi,bj,k, dTtracerLev(k),
351         I           vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,
352         O           af, myThid )
353             ENDIF
354            ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN
355              CALL GAD_OS7MP_ADV_Y( bi,bj,k, dTtracerLev(k),
356         I            vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,
357         O            af, myThid )
358            ELSE
359              STOP 'GAD_CALC_RHS: Bad advectionScheme (Y)'
360            ENDIF
361            DO j=1-Oly,sNy+Oly
362             DO i=1-Olx,sNx+Olx
363              fMer(i,j) = fMer(i,j) + af(i,j)
364             ENDDO
365            ENDDO
366    #ifdef ALLOW_DIAGNOSTICS
367            IF ( useDiagnostics ) THEN
368              diagName = 'ADVy'//diagSufx
369              CALL DIAGNOSTICS_FILL(af,diagName, k,1, 2,bi,bj, myThid)
370            ENDIF
371    #endif
372        ENDIF        ENDIF
       DO j=1-Oly,sNy+Oly  
        DO i=1-Olx,sNx+Olx  
         fMer(i,j) = fMer(i,j) + af(i,j)  
        ENDDO  
       ENDDO  
373    
374  C-    Diffusive flux in Y  C-    Diffusive flux in Y
375        IF (diffKh.NE.0.) THEN        IF (diffKh.NE.0.) THEN
# Line 194  C-    Diffusive flux in Y Line 377  C-    Diffusive flux in Y
377        ELSE        ELSE
378         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
379          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
380           df(i,j) = 0.           df(i,j) = 0. _d 0
381          ENDDO          ENDDO
382         ENDDO         ENDDO
383        ENDIF        ENDIF
384    
385    C-    Add bi-harmonic flux in Y
386          IF (diffK4 .NE. 0.) THEN
387           CALL GAD_BIHARM_Y(bi,bj,k,yA,df4,diffK4,df,myThid)
388          ENDIF
389    
390  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
391  C-    GM/Redi flux in Y  C-    GM/Redi flux in Y
392        IF (useGMRedi) THEN        IF (useGMRedi) THEN
393         CALL GMREDI_YTRANSPORT(  C *note* should update GMREDI_YTRANSPORT to set df  *aja*
394  C *note* should update GMREDI_YTRANSPORT to use localT and set df  *aja*          IF ( applyAB_onTracer ) THEN
395       I     iMin,iMax,jMin,jMax,bi,bj,K,            CALL GMREDI_YTRANSPORT(
396       I     yA,Tracer,       I         iMin,iMax,jMin,jMax,bi,bj,k,
397       U     df,       I         yA,TracerN,tracerIdentity,
398       I     myThid)       U         df,
399         I         myThid)
400            ELSE
401              CALL GMREDI_YTRANSPORT(
402         I         iMin,iMax,jMin,jMax,bi,bj,k,
403         I         yA,TracAB, tracerIdentity,
404         U         df,
405         I         myThid)
406            ENDIF
407        ENDIF        ENDIF
408  #endif  #endif
409    C     anelastic: advect.fluxes are scaled by rhoFac but hor.diff. flx are not
410        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
411         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
412          fMer(i,j) = fMer(i,j) + df(i,j)          fMer(i,j) = fMer(i,j) + df(i,j)*rhoFacC(k)
413         ENDDO         ENDDO
414        ENDDO        ENDDO
415    
416  C-    Bi-harmonic flux in Y  #ifdef ALLOW_DIAGNOSTICS
417        IF (diffK4 .NE. 0.) THEN  C-    Diagnostics of Tracer flux in Y dir (mainly Diffusive terms),
418         CALL GAD_BIHARM_Y(bi,bj,k,yA,df4,diffK4,df,myThid)  C       excluding advective terms:
419         DO j=1-Oly,sNy+Oly        IF ( useDiagnostics .AND.
420          DO i=1-Olx,sNx+Olx       &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. useGMRedi) ) THEN
421           fMer(i,j) = fMer(i,j) + df(i,j)            diagName = 'DFyE'//diagSufx
422          ENDDO            CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
        ENDDO  
423        ENDIF        ENDIF
424    #endif
425    
426  C--   Initialize net flux in R  C--   Compute vertical flux fVerT(kUp) at interface k (between k-1 & k):
       DO j=1-Oly,sNy+Oly  
        DO i=1-Olx,sNx+Olx  
         fVerT(i,j,kUp) = 0.  
        ENDDO  
       ENDDO  
   
427  C-    Advective flux in R  C-    Advective flux in R
428  C     Note: wVel needs to be masked  #ifdef ALLOW_AIM
429        IF (K.GE.2) THEN  C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr
430          IF (calcAdvection .AND. .NOT.implicitAdvection .AND. k.GE.2 .AND.
431         &     (.NOT.useAIM .OR.tracerIdentity.NE.GAD_SALINITY .OR.k.LT.Nr)
432         &   ) THEN
433    #else
434          IF (calcAdvection .AND. .NOT.implicitAdvection .AND. k.GE.2) THEN
435    #endif
436  C-    Compute vertical advective flux in the interior:  C-    Compute vertical advective flux in the interior:
437         IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN          IF (vertAdvecScheme.EQ.ENUM_CENTERED_2ND) THEN
438          CALL GAD_C2_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)            CALL GAD_C2_ADV_R(bi,bj,k,rTrans,TracAB,af,myThid)
439         ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF ( vertAdvecScheme.EQ.ENUM_UPWIND_1RST
440          CALL GAD_FLUXLIMIT_ADV_R(       &          .OR. vertAdvecScheme.EQ.ENUM_DST2 ) THEN
441       &       bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)            CALL GAD_DST2U1_ADV_R( bi,bj,k, vertAdvecScheme,
442         ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN       I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),
443          CALL GAD_U3_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)       O         af, myThid )
444         ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN          ELSEIF (vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN
445          CALL GAD_C4_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)            CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k,
446         ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN       I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),
447  c       CALL GAD_DST3_ADV_R(       O         af, myThid )
448  c    &       bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)          ELSEIF (vertAdvecScheme.EQ.ENUM_UPWIND_3RD ) THEN
449          STOP 'GAD_CALC_RHS: GAD_DST3_ADV_R not coded yet'            CALL GAD_U3_ADV_R(bi,bj,k,rTrans,TracAB,af,myThid)
450         ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN          ELSEIF (vertAdvecScheme.EQ.ENUM_CENTERED_4TH) THEN
451  c       CALL GAD_DST3FL_ADV_R(            CALL GAD_C4_ADV_R(bi,bj,k,rTrans,TracAB,af,myThid)
452  c    &       bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)          ELSEIF (vertAdvecScheme.EQ.ENUM_DST3 ) THEN
453          STOP 'GAD_CALC_RHS: GAD_DST3FL_ADV_R not coded yet'            CALL GAD_DST3_ADV_R( bi,bj,k,
454         ELSE       I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),
455          STOP 'GAD_CALC_RHS: Bad advectionScheme (R)'       O         af, myThid )
456         ENDIF          ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
457  C-    Surface "correction" term at k>1 :  cph This block is to trick the adjoint:
458         DO j=1-Oly,sNy+Oly  cph IF inAdExact=.FALSE., we want to use DST3
459          DO i=1-Olx,sNx+Olx  cph with limiters in forward, but without limiters in reverse.
460           af(i,j) = af(i,j)            IF ( inAdMode ) THEN
461       &           + (maskC(i,j,k,bi,bj)-maskC(i,j,k-1,bi,bj))*             CALL GAD_DST3_ADV_R( bi,bj,k,
462       &             rTrans(i,j)*Tracer(i,j,k,bi,bj)       I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),
463         O         af, myThid )
464              ELSE
465               CALL GAD_DST3FL_ADV_R( bi,bj,k,
466         I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),
467         O         af, myThid )
468              ENDIF
469            ELSEIF (vertAdvecScheme.EQ.ENUM_OS7MP ) THEN
470               CALL GAD_OS7MP_ADV_R( bi,bj,k,
471         I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),
472         O         af, myThid )
473            ELSE
474              STOP 'GAD_CALC_RHS: Bad vertAdvecScheme (R)'
475            ENDIF
476    C-     add the advective flux to fVerT
477            DO j=1-Oly,sNy+Oly
478             DO i=1-Olx,sNx+Olx
479              fVerT(i,j,kUp) = fVerT(i,j,kUp) + af(i,j)
480             ENDDO
481          ENDDO          ENDDO
482         ENDDO  #ifdef ALLOW_DIAGNOSTICS
483        ELSE          IF ( useDiagnostics ) THEN
484  C-    Surface "correction" term at k=1 :            diagName = 'ADVr'//diagSufx
485         DO j=1-Oly,sNy+Oly            CALL DIAGNOSTICS_FILL(af,diagName, k,1, 2,bi,bj, myThid)
486          DO i=1-Olx,sNx+Olx  C- note: needs to explicitly increment the counter since DIAGNOSTICS_FILL
487           af(i,j) = rTrans(i,j)*Tracer(i,j,k,bi,bj)  C        does it only if k=1 (never the case here)
488          ENDDO            IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT(diagName,bi,bj,myThid)
489         ENDDO          ENDIF
490    #endif
491        ENDIF        ENDIF
 C-    add the advective flux to fVerT  
       DO j=1-Oly,sNy+Oly  
        DO i=1-Olx,sNx+Olx  
         fVerT(i,j,kUp) = fVerT(i,j,kUp) + afFacT*af(i,j)  
        ENDDO  
       ENDDO  
492    
493  C-    Diffusive flux in R  C-    Diffusive flux in R
494  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
# Line 286  C           boundary condition. Line 496  C           boundary condition.
496        IF (implicitDiffusion) THEN        IF (implicitDiffusion) THEN
497         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
498          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
499           df(i,j) = 0.           df(i,j) = 0. _d 0
500          ENDDO          ENDDO
501         ENDDO         ENDDO
502        ELSE        ELSE
503         CALL GAD_DIFF_R(bi,bj,k,KappaRT,tracer,df,myThid)         IF ( applyAB_onTracer ) THEN
504             CALL GAD_DIFF_R(bi,bj,k,KappaR,TracerN,df,myThid)
505           ELSE
506             CALL GAD_DIFF_R(bi,bj,k,KappaR,TracAB, df,myThid)
507           ENDIF
508        ENDIF        ENDIF
 c     DO j=1-Oly,sNy+Oly  
 c      DO i=1-Olx,sNx+Olx  
 c       fVerT(i,j,kUp) = fVerT(i,j,kUp) + dfFacT*df(i,j)*maskUp(i,j)  
 c      ENDDO  
 c     ENDDO  
509    
510  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
511  C-    GM/Redi flux in R  C-    GM/Redi flux in R
512        IF (useGMRedi) THEN        IF (useGMRedi) THEN
513  C *note* should update GMREDI_RTRANSPORT to set df  *aja*  C *note* should update GMREDI_RTRANSPORT to set df  *aja*
514         CALL GMREDI_RTRANSPORT(          IF ( applyAB_onTracer ) THEN
515       I     iMin,iMax,jMin,jMax,bi,bj,K,            CALL GMREDI_RTRANSPORT(
516       I     Tracer,       I         iMin,iMax,jMin,jMax,bi,bj,k,
517       U     df,       I         TracerN,tracerIdentity,
518       I     myThid)       U         df,
519  c      DO j=1-Oly,sNy+Oly       I         myThid)
520  c       DO i=1-Olx,sNx+Olx          ELSE
521  c        fVerT(i,j,kUp) = fVerT(i,j,kUp) + dfFacT*df(i,j)*maskUp(i,j)            CALL GMREDI_RTRANSPORT(
522  c       ENDDO       I         iMin,iMax,jMin,jMax,bi,bj,k,
523  c      ENDDO       I         TracAB, tracerIdentity,
524         U         df,
525         I         myThid)
526            ENDIF
527        ENDIF        ENDIF
528  #endif  #endif
529    
530        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
531         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
532          fVerT(i,j,kUp) = fVerT(i,j,kUp) + dfFacT*df(i,j)*maskUp(i,j)          fVerT(i,j,kUp) = fVerT(i,j,kUp) + df(i,j)*maskUp(i,j)
533         ENDDO         ENDDO
534        ENDDO        ENDDO
535    
536    #ifdef ALLOW_DIAGNOSTICS
537    C-    Diagnostics of Tracer flux in R dir (mainly Diffusive terms),
538    C       Explicit terms only & excluding advective terms:
539          IF ( useDiagnostics .AND.
540         &    (.NOT.implicitDiffusion .OR. useGMRedi) ) THEN
541              diagName = 'DFrE'//diagSufx
542              CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
543          ENDIF
544    #endif
545    
546  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
547  C-    Add non local KPP transport term (ghat) to diffusive T flux.  C-    Set non local KPP transport term (ghat):
548        IF (useKPP) THEN        IF ( useKPP .AND. k.GE.2 ) THEN
549         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
550          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
551           df(i,j) = 0.           df(i,j) = 0. _d 0
552          ENDDO          ENDDO
553         ENDDO         ENDDO
554         IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN         IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
 C *note* should update KPP_TRANSPORT_T to set df  *aja*  
555          CALL KPP_TRANSPORT_T(          CALL KPP_TRANSPORT_T(
556       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,
557       I     KappaRT,       O     df )
      U     df )  
558         ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN         ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
559          CALL KPP_TRANSPORT_S(          CALL KPP_TRANSPORT_S(
560       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,
561       I     KappaRT,       O     df )
562       U     df )  #ifdef ALLOW_PTRACERS
563           ELSEIF (tracerIdentity .GE. GAD_TR1) THEN
564            CALL KPP_TRANSPORT_PTR(
565         I     iMin,iMax,jMin,jMax,bi,bj,k,km1,
566         I     tracerIdentity-GAD_TR1+1,
567         O     df )
568    #endif
569         ELSE         ELSE
570            PRINT*,'invalid tracer indentity: ', tracerIdentity
571          STOP 'GAD_CALC_RHS: Ooops'          STOP 'GAD_CALC_RHS: Ooops'
572         ENDIF         ENDIF
573         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
574          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
575           fVerT(i,j,kUp) = fVerT(i,j,kUp) + dfFacT*df(i,j)*maskUp(i,j)           fVerT(i,j,kUp) = fVerT(i,j,kUp)
576         &                  + df(i,j)*maskUp(i,j)*rhoFacF(k)
577          ENDDO          ENDDO
578         ENDDO         ENDDO
579        ENDIF        ENDIF
580  #endif  #endif
581    
582  C--   Divergence of fluxes  C--   Divergence of fluxes
583        DO j=1-Oly,sNy+Oly  C     Anelastic: scale vertical fluxes by rhoFac and leave Horizontal fluxes unchanged
584         DO i=1-Olx,sNx+Olx        DO j=1-Oly,sNy+Oly-1
585          gTracer(i,j,k,bi,bj)=         DO i=1-Olx,sNx+Olx-1
586            gTracer(i,j,k,bi,bj)=gTracer(i,j,k,bi,bj)
587       &   -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &   -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
588       &    *recip_rA(i,j,bi,bj)       &   *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)*recip_rhoFacC(k)
589       &    *(       &   *( (fZon(i+1,j)-fZon(i,j))
590       &    +( fZon(i+1,j)-fZon(i,j) )       &     +(fMer(i,j+1)-fMer(i,j))
591       &    +( fMer(i,j+1)-fMer(i,j) )       &     +(fVerT(i,j,kDown)-fVerT(i,j,kUp))*rkSign
592       &    +( fVerT(i,j,kUp)-fVerT(i,j,kDown) )*rkFac       &     -localT(i,j)*( (uTrans(i+1,j)-uTrans(i,j))
593         &                   +(vTrans(i,j+1)-vTrans(i,j))
594         &                   +(rTransKp1(i,j)-rTrans(i,j))*rAdvFac
595         &                  )*advFac
596       &    )       &    )
597         ENDDO         ENDDO
598        ENDDO        ENDDO
599    
600    #ifdef ALLOW_DEBUG
601          IF ( debugLevel .GE. debLevB
602         &   .AND. tracerIdentity.EQ.GAD_TEMPERATURE
603         &   .AND. k.EQ.2 .AND. myIter.EQ.1+nIter0
604         &   .AND. nPx.EQ.1 .AND. nPy.EQ.1
605         &   .AND. useCubedSphereExchange ) THEN
606            CALL DEBUG_CS_CORNER_UV( ' fZon,fMer from GAD_CALC_RHS',
607         &             fZon,fMer, k, standardMessageUnit,bi,bj,myThid )
608          ENDIF
609    #endif /* ALLOW_DEBUG */
610    
611        RETURN        RETURN
612        END        END
Legend:
Removed from v.1.6  
changed lines
  Added in v.1.44
  ViewVC Help
Powered by ViewVC 1.1.22