/[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.4 by adcroft, Tue Sep 4 14:53:11 2001 UTC revision 1.34 by jmc, Thu Feb 17 00:03:07 2005 UTC
# Line 3  C $Name$ Line 3  C $Name$
3    
4  #include "GAD_OPTIONS.h"  #include "GAD_OPTIONS.h"
5    
6    CBOP
7    C !ROUTINE: GAD_CALC_RHS
8    
9    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,uTrans,vTrans,rTrans,rTransKp1,maskUp,
13       I           diffKh, diffK4, KappaRT, Tracer,       I           uVel, vVel, wVel,
14       I           tracerIdentity, advectionScheme,       I           diffKh, diffK4, KappaR, Tracer,
15         I           tracerIdentity, advectionScheme, vertAdvecScheme,
16         I           calcAdvection, implicitAdvection,
17       U           fVerT, gTracer,       U           fVerT, gTracer,
18       I           myThid )       I           myTime, myIter, myThid )
19  C     /==========================================================\  
20  C     | SUBROUTINE GAD_CALC_RHS                                  |  C !DESCRIPTION:
21  C     |==========================================================|  C Calculates the tendancy of a tracer due to advection and diffusion.
22  C     \==========================================================/  C It calculates the fluxes in each direction indepentently and then
23        IMPLICIT NONE  C sets the tendancy 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 tendancy 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 tendancy is assumed to contain data on entry.
38    
39  C     == GLobal variables ==  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 kup              :: index into 2 1/2D array, toggles between 1|2
58    C kdown            :: index into 2 1/2D array, toggles between 2|1
59    C kp1              :: =k+1 for k<Nr, =Nr for k=Nr
60    C xA,yA            :: areas of X and Y face of tracer cells
61    C uTrans,vTrans    :: 2-D arrays of volume transports at U,V points
62    C rTrans           :: 2-D arrays of volume transports at W points
63    C rTransKp1        :: 2-D array of volume trans at W pts, interf k+1
64    C maskUp           :: 2-D array for mask at W points
65    C uVel,vVel,wVel   :: 3 components of the velcity field (3-D array)
66    C diffKh           :: horizontal diffusion coefficient
67    C diffK4           :: bi-harmonic diffusion coefficient
68    C KappaR           :: 2-D array for vertical diffusion coefficient, interf k
69    C Tracer           :: tracer field
70    C tracerIdentity   :: tracer identifier (required for KPP,GM)
71    C advectionScheme  :: advection scheme to use (Horizontal plane)
72    C vertAdvecScheme  :: advection scheme to use (Vertical direction)
73    C calcAdvection    :: =False if Advec computed with multiDim scheme
74    C implicitAdvection:: =True if vertical Advec computed implicitly
75    C myTime           :: current time
76    C myIter           :: iteration number
77    C myThid           :: thread number
78        INTEGER bi,bj,iMin,iMax,jMin,jMax        INTEGER bi,bj,iMin,iMax,jMin,jMax
79          INTEGER k,kUp,kDown,kM1
80        _RS xA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS xA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
81        _RS yA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS yA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
82        _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
83        _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
84        _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
85          _RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
86        _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
87          _RL uVel  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
88          _RL vVel  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
89          _RL wVel  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
90        _RL diffKh, diffK4        _RL diffKh, diffK4
91        _RL KappaRT(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL KappaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
92        _RL Tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL Tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
93        INTEGER tracerIdentity        INTEGER tracerIdentity
94        INTEGER advectionScheme        INTEGER advectionScheme, vertAdvecScheme
95        _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        LOGICAL calcAdvection
96          LOGICAL implicitAdvection
97          _RL     myTime
98          INTEGER myIter, myThid
99    
100    C !OUTPUT PARAMETERS: ==================================================
101    C gTracer          :: tendancy array
102    C fVerT            :: 2 1/2D arrays for vertical advective flux
103        _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)
104        INTEGER myThid        _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
105    
106  C     == Local variables ==  C !LOCAL VARIABLES: ====================================================
107  C     I, J, K - Loop counters  C i,j              :: loop indices
108    C df4              :: used for storing del^2 T for bi-harmonic term
109    C fZon             :: zonal flux
110    C fMer             :: meridional flux
111    C af               :: advective flux
112    C df               :: diffusive flux
113    C localT           :: local copy of tracer field
114    #ifdef ALLOW_DIAGNOSTICS
115          CHARACTER*8 diagName
116          CHARACTER*4 GAD_DIAG_SUFX, diagSufx
117          EXTERNAL    GAD_DIAG_SUFX
118    #endif
119        INTEGER i,j        INTEGER i,j
       LOGICAL TOP_LAYER  
       _RL afFacT, dfFacT  
120        _RL df4   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL df4   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
121        _RL fZon  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fZon  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
122        _RL fMer  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fMer  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
123        _RL af    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL af    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
124        _RL df    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL df    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
125        _RL localT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL localT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
126          _RL advFac, rAdvFac
127    CEOP
128    
129  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
130  C--   only the kUp part of fverT is set in this subroutine  C--   only the kUp part of fverT is set in this subroutine
131  C--   the kDown is still required  C--   the kDown is still required
132        fVerT(1,1,kDown) = fVerT(1,1,kDown)        fVerT(1,1,kDown) = fVerT(1,1,kDown)
133  #endif  #endif
134    
135    #ifdef ALLOW_DIAGNOSTICS
136    C--   Set diagnostic suffix for the current tracer
137          IF ( useDiagnostics ) THEN
138            diagSufx = GAD_DIAG_SUFX( tracerIdentity, myThid )
139          ENDIF
140    #endif
141    
142          advFac  = 0. _d 0
143          IF (calcAdvection) advFac = 1. _d 0
144          rAdvFac = rkFac*advFac
145          IF (implicitAdvection) rAdvFac = 0. _d 0
146    
147        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
148         DO i=1-OLx,sNx+OLx         DO i=1-OLx,sNx+OLx
149          fZon(i,j)      = 0.0          fZon(i,j)      = 0. _d 0
150          fMer(i,j)      = 0.0          fMer(i,j)      = 0. _d 0
151          fVerT(i,j,kUp) = 0.0          fVerT(i,j,kUp) = 0. _d 0
152            df(i,j)        = 0. _d 0
153            df4(i,j)       = 0. _d 0
154         ENDDO         ENDDO
155        ENDDO        ENDDO
156    
       afFacT = 1. _d 0  
       dfFacT = 1. _d 0  
       TOP_LAYER = K .EQ. 1  
   
157  C--   Make local copy of tracer array  C--   Make local copy of tracer array
158        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
159         DO i=1-OLx,sNx+OLx         DO i=1-OLx,sNx+OLx
# Line 78  C--   Make local copy of tracer array Line 161  C--   Make local copy of tracer array
161         ENDDO         ENDDO
162        ENDDO        ENDDO
163    
164    C--   Unless we have already calculated the advection terms we initialize
165    C     the tendency to zero.
166    C     <== now done earlier at the beginning of thermodynamics.
167    c     IF (calcAdvection) THEN
168    c      DO j=1-Oly,sNy+Oly
169    c       DO i=1-Olx,sNx+Olx
170    c        gTracer(i,j,k,bi,bj)=0. _d 0
171    c       ENDDO
172    c      ENDDO
173    c     ENDIF
174    
175  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
176        IF (diffK4 .NE. 0.) THEN        IF (diffK4 .NE. 0.) THEN
# Line 89  C--   Pre-calculate del^2 T if bi-harmon Line 182  C--   Pre-calculate del^2 T if bi-harmon
182  C--   Initialize net flux in X direction  C--   Initialize net flux in X direction
183        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
184         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
185          fZon(i,j) = 0.          fZon(i,j) = 0. _d 0
186         ENDDO         ENDDO
187        ENDDO        ENDDO
188    
189  C-    Advective flux in X  C-    Advective flux in X
190        IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN        IF (calcAdvection) THEN
191         CALL GAD_C2_ADV_X(bi,bj,k,uTrans,localT,af,myThid)          IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
192        ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN            CALL GAD_C2_ADV_X(bi,bj,k,uTrans,localT,af,myThid)
193         CALL GAD_FLUXLIMIT_ADV_X(          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
194       &      bi,bj,k,deltaTtracer,uTrans,uVel,localT,af,myThid)            CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, dTtracerLev(k),
195        ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN       I            uTrans, uVel, maskW(1-Olx,1-Oly,k,bi,bj), localT,
196         CALL GAD_U3_ADV_X(bi,bj,k,uTrans,localT,af,myThid)       O            af, myThid )
197        ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN          ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
198         CALL GAD_C4_ADV_X(bi,bj,k,uTrans,localT,af,myThid)            CALL GAD_U3_ADV_X(bi,bj,k,uTrans,localT,af,myThid)
199        ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN          ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
200         CALL GAD_DST3_ADV_X(            CALL GAD_C4_ADV_X(bi,bj,k,uTrans,localT,af,myThid)
201       &       bi,bj,k,deltaTtracer,uTrans,uVel,localT,af,myThid)          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
202        ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN            CALL GAD_DST3_ADV_X( bi,bj,k, dTtracerLev(k),
203         CALL GAD_DST3FL_ADV_X(       I            uTrans, uVel, maskW(1-Olx,1-Oly,k,bi,bj), localT,
204       &       bi,bj,k,deltaTtracer,uTrans,uVel,localT,af,myThid)       O            af, myThid )
205        ELSE          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
206         STOP 'GAD_CALC_RHS: Bad advectionScheme (X)'            CALL GAD_DST3FL_ADV_X( bi,bj,k, dTtracerLev(k),
207         I            uTrans, uVel, maskW(1-Olx,1-Oly,k,bi,bj), localT,
208         O            af, myThid )
209            ELSE
210             STOP 'GAD_CALC_RHS: Bad advectionScheme (X)'
211            ENDIF
212            DO j=1-Oly,sNy+Oly
213             DO i=1-Olx,sNx+Olx
214              fZon(i,j) = fZon(i,j) + af(i,j)
215             ENDDO
216            ENDDO
217    #ifdef ALLOW_DIAGNOSTICS
218            IF ( useDiagnostics ) THEN
219              diagName = 'ADVx'//diagSufx
220              CALL DIAGNOSTICS_FILL(af,diagName, k,1, 2,bi,bj, myThid)
221            ENDIF
222    #endif
223        ENDIF        ENDIF
       DO j=jMin,jMax  
        DO i=iMin,iMax  
         fZon(i,j) = fZon(i,j) + af(i,j)  
        ENDDO  
       ENDDO  
224    
225  C-    Diffusive flux in X  C-    Diffusive flux in X
226        IF (diffKh.NE.0.) THEN        IF (diffKh.NE.0.) THEN
227         CALL GAD_DIFF_X(bi,bj,k,xA,diffKh,localT,df,myThid)         CALL GAD_DIFF_X(bi,bj,k,xA,diffKh,localT,df,myThid)
228        ELSE        ELSE
229         DO j=jMin,jMax         DO j=1-Oly,sNy+Oly
230          DO i=iMin,iMax          DO i=1-Olx,sNx+Olx
231           df(i,j) = 0.           df(i,j) = 0. _d 0
232          ENDDO          ENDDO
233         ENDDO         ENDDO
234        ENDIF        ENDIF
235    
236    C-    Add bi-harmonic diffusive flux in X
237          IF (diffK4 .NE. 0.) THEN
238           CALL GAD_BIHARM_X(bi,bj,k,xA,df4,diffK4,df,myThid)
239          ENDIF
240    
241  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
242  C-    GM/Redi flux in X  C-    GM/Redi flux in X
243        IF (useGMRedi) THEN        IF (useGMRedi) THEN
244  C *note* should update GMREDI_XTRANSPORT to use localT and set df  *aja*  C *note* should update GMREDI_XTRANSPORT to use localT and set df  *aja*
245          CALL GMREDI_XTRANSPORT(          CALL GMREDI_XTRANSPORT(
246       I     iMin,iMax,jMin,jMax,bi,bj,K,       I     iMin,iMax,jMin,jMax,bi,bj,K,
247       I     xA,Tracer,       I     xA,Tracer,tracerIdentity,
248       U     df,       U     df,
249       I     myThid)       I     myThid)
250        ENDIF        ENDIF
251  #endif  #endif
252        DO j=jMin,jMax        DO j=1-Oly,sNy+Oly
253         DO i=iMin,iMax         DO i=1-Olx,sNx+Olx
254          fZon(i,j) = fZon(i,j) + df(i,j)          fZon(i,j) = fZon(i,j) + df(i,j)
255         ENDDO         ENDDO
256        ENDDO        ENDDO
257    
258  C-    Bi-harmonic duffusive flux in X  #ifdef ALLOW_DIAGNOSTICS
259        IF (diffK4 .NE. 0.) THEN  C-    Diagnostics of Tracer flux in X dir (mainly Diffusive term),
260         CALL GAD_BIHARM_X(bi,bj,k,xA,df4,diffK4,df,myThid)  C       excluding advective terms:
261         DO j=jMin,jMax        IF ( useDiagnostics .AND.
262          DO i=iMin,iMax       &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. useGMRedi) ) THEN
263           fZon(i,j) = fZon(i,j) + df(i,j)            diagName = 'DIFx'//diagSufx
264          ENDDO            CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
        ENDDO  
265        ENDIF        ENDIF
266    #endif
267    
268  C--   Initialize net flux in Y direction  C--   Initialize net flux in Y direction
269        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
270         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
271          fMer(i,j) = 0.          fMer(i,j) = 0. _d 0
272         ENDDO         ENDDO
273        ENDDO        ENDDO
274    
275  C-    Advective flux in Y  C-    Advective flux in Y
276        IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN        IF (calcAdvection) THEN
277         CALL GAD_C2_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)          IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
278        ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN            CALL GAD_C2_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)
279         CALL GAD_FLUXLIMIT_ADV_Y(          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
280       &       bi,bj,k,deltaTtracer,vTrans,vVel,localT,af,myThid)            CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, dTtracerLev(k),
281        ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN       I            vTrans, vVel, maskS(1-Olx,1-Oly,k,bi,bj), localT,
282         CALL GAD_U3_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)       O            af, myThid )
283        ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN          ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
284         CALL GAD_C4_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)            CALL GAD_U3_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)
285        ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN          ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
286         CALL GAD_DST3_ADV_Y(            CALL GAD_C4_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)
287       &       bi,bj,k,deltaTtracer,vTrans,vVel,localT,af,myThid)          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
288        ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN            CALL GAD_DST3_ADV_Y( bi,bj,k, dTtracerLev(k),
289         CALL GAD_DST3FL_ADV_Y(       I            vTrans, vVel, maskS(1-Olx,1-Oly,k,bi,bj), localT,
290       &       bi,bj,k,deltaTtracer,vTrans,vVel,localT,af,myThid)       O            af, myThid )
291        ELSE          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
292         STOP 'GAD_CALC_RHS: Bad advectionScheme (Y)'            CALL GAD_DST3FL_ADV_Y( bi,bj,k, dTtracerLev(k),
293         I            vTrans, vVel, maskS(1-Olx,1-Oly,k,bi,bj), localT,
294         O            af, myThid )
295            ELSE
296              STOP 'GAD_CALC_RHS: Bad advectionScheme (Y)'
297            ENDIF
298            DO j=1-Oly,sNy+Oly
299             DO i=1-Olx,sNx+Olx
300              fMer(i,j) = fMer(i,j) + af(i,j)
301             ENDDO
302            ENDDO
303    #ifdef ALLOW_DIAGNOSTICS
304            IF ( useDiagnostics ) THEN
305              diagName = 'ADVy'//diagSufx
306              CALL DIAGNOSTICS_FILL(af,diagName, k,1, 2,bi,bj, myThid)
307            ENDIF
308    #endif
309        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  
310    
311  C-    Diffusive flux in Y  C-    Diffusive flux in Y
312        IF (diffKh.NE.0.) THEN        IF (diffKh.NE.0.) THEN
# Line 194  C-    Diffusive flux in Y Line 314  C-    Diffusive flux in Y
314        ELSE        ELSE
315         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
316          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
317           df(i,j) = 0.           df(i,j) = 0. _d 0
318          ENDDO          ENDDO
319         ENDDO         ENDDO
320        ENDIF        ENDIF
321    
322    C-    Add bi-harmonic flux in Y
323          IF (diffK4 .NE. 0.) THEN
324           CALL GAD_BIHARM_Y(bi,bj,k,yA,df4,diffK4,df,myThid)
325          ENDIF
326    
327  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
328  C-    GM/Redi flux in Y  C-    GM/Redi flux in Y
329        IF (useGMRedi) THEN        IF (useGMRedi) THEN
        CALL GMREDI_YTRANSPORT(  
330  C *note* should update GMREDI_YTRANSPORT to use localT and set df  *aja*  C *note* should update GMREDI_YTRANSPORT to use localT and set df  *aja*
331           CALL GMREDI_YTRANSPORT(
332       I     iMin,iMax,jMin,jMax,bi,bj,K,       I     iMin,iMax,jMin,jMax,bi,bj,K,
333       I     yA,Tracer,       I     yA,Tracer,tracerIdentity,
334       U     df,       U     df,
335       I     myThid)       I     myThid)
336        ENDIF        ENDIF
# Line 216  C *note* should update GMREDI_YTRANSPORT Line 341  C *note* should update GMREDI_YTRANSPORT
341         ENDDO         ENDDO
342        ENDDO        ENDDO
343    
344  C-    Bi-harmonic flux in Y  #ifdef ALLOW_DIAGNOSTICS
345        IF (diffK4 .NE. 0.) THEN  C-    Diagnostics of Tracer flux in Y dir (mainly Diffusive terms),
346         CALL GAD_BIHARM_Y(bi,bj,k,yA,df4,diffK4,df,myThid)  C       excluding advective terms:
347         DO j=1-Oly,sNy+Oly        IF ( useDiagnostics .AND.
348          DO i=1-Olx,sNx+Olx       &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. useGMRedi) ) THEN
349           fMer(i,j) = fMer(i,j) + df(i,j)            diagName = 'DIFy'//diagSufx
350          ENDDO            CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
        ENDDO  
351        ENDIF        ENDIF
352    #endif
353    
354  C--   Initialize net flux in R  C--   Compute vertical flux fVerT(kUp) at interface k (between k-1 & k):
       DO j=jMin,jMax  
        DO i=iMin,iMax  
         fVerT(i,j,kUp) = 0.  
        ENDDO  
       ENDDO  
   
355  C-    Advective flux in R  C-    Advective flux in R
356  C     Note: wVel needs to be masked  #ifdef ALLOW_AIM
357        IF (K.GE.2) THEN  C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr
358          IF (calcAdvection .AND. .NOT.implicitAdvection .AND. K.GE.2 .AND.
359         &     (.NOT.useAIM .OR.tracerIdentity.NE.GAD_SALINITY .OR.K.LT.Nr)
360         &   ) THEN
361    #else
362          IF (calcAdvection .AND. .NOT.implicitAdvection .AND. K.GE.2) THEN
363    #endif
364  C-    Compute vertical advective flux in the interior:  C-    Compute vertical advective flux in the interior:
365         IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN          IF (vertAdvecScheme.EQ.ENUM_CENTERED_2ND) THEN
366          CALL GAD_C2_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)            CALL GAD_C2_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)
367         ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF (vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN
368          CALL GAD_FLUXLIMIT_ADV_R(            CALL GAD_FLUXLIMIT_ADV_R(
369       &       bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)       &         bi,bj,k,dTtracerLev(k),rTrans,wVel,tracer,af,myThid)
370         ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN          ELSEIF (vertAdvecScheme.EQ.ENUM_UPWIND_3RD ) THEN
371          CALL GAD_U3_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)            CALL GAD_U3_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)
372         ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN          ELSEIF (vertAdvecScheme.EQ.ENUM_CENTERED_4TH) THEN
373          CALL GAD_C4_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)            CALL GAD_C4_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)
374         ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN          ELSEIF (vertAdvecScheme.EQ.ENUM_DST3 ) THEN
375  c       CALL GAD_DST3_ADV_R(            CALL GAD_DST3_ADV_R(
376  c    &       bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)       &         bi,bj,k,dTtracerLev(k),rTrans,wVel,tracer,af,myThid)
377          STOP 'GAD_CALC_RHS: GAD_DST3_ADV_R not coded yet'          ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
378         ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN            CALL GAD_DST3FL_ADV_R(
379  c       CALL GAD_DST3FL_ADV_R(       &         bi,bj,k,dTtracerLev(k),rTrans,wVel,tracer,af,myThid)
380  c    &       bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)          ELSE
381          STOP 'GAD_CALC_RHS: GAD_DST3FL_ADV_R not coded yet'            STOP 'GAD_CALC_RHS: Bad vertAdvecScheme (R)'
382         ELSE          ENDIF
383          STOP 'GAD_CALC_RHS: Bad advectionScheme (R)'  C-     add the advective flux to fVerT
384         ENDIF          DO j=1-Oly,sNy+Oly
385  C-    Surface "correction" term at k>1 :           DO i=1-Olx,sNx+Olx
386         DO j=1-Oly,sNy+Oly            fVerT(i,j,kUp) = fVerT(i,j,kUp) + af(i,j)
387          DO i=1-Olx,sNx+Olx           ENDDO
          af(i,j) = af(i,j)  
      &           + (maskC(i,j,k,bi,bj)-maskC(i,j,k-1,bi,bj))*  
      &             rTrans(i,j)*Tracer(i,j,k,bi,bj)  
         ENDDO  
        ENDDO  
       ELSE  
 C-    Surface "correction" term at k=1 :  
        DO j=1-Oly,sNy+Oly  
         DO i=1-Olx,sNx+Olx  
          af(i,j) = rTrans(i,j)*Tracer(i,j,k,bi,bj)  
388          ENDDO          ENDDO
389         ENDDO  #ifdef ALLOW_DIAGNOSTICS
390            IF ( useDiagnostics ) THEN
391              diagName = 'ADVr'//diagSufx
392              CALL DIAGNOSTICS_FILL(af,diagName, k,1, 2,bi,bj, myThid)
393    C- note: needs to explicitly increment the counter since DIAGNOSTICS_FILL
394    C        does it only if k=1 (never the case here)
395              IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT(diagName,bi,bj,myThid)
396            ENDIF
397    #endif
398        ENDIF        ENDIF
 C-    add the advective flux to fVerT  
       DO j=jMin,jMax  
        DO i=iMin,iMax  
         fVerT(i,j,kUp) = fVerT(i,j,kUp) + afFacT*af(i,j)  
        ENDDO  
       ENDDO  
399    
400  C-    Diffusive flux in R  C-    Diffusive flux in R
401  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
402  C           boundary condition.  C           boundary condition.
403        IF (implicitDiffusion) THEN        IF (implicitDiffusion) THEN
404         DO j=jMin,jMax         DO j=1-Oly,sNy+Oly
405          DO i=iMin,iMax          DO i=1-Olx,sNx+Olx
406           df(i,j) = 0.           df(i,j) = 0. _d 0
407          ENDDO          ENDDO
408         ENDDO         ENDDO
409        ELSE        ELSE
410         CALL GAD_DIFF_R(bi,bj,k,KappaRT,tracer,df,myThid)         CALL GAD_DIFF_R(bi,bj,k,KappaR,tracer,df,myThid)
411        ENDIF        ENDIF
 c     DO j=jMin,jMax  
 c      DO i=iMin,iMax  
 c       fVerT(i,j,kUp) = fVerT(i,j,kUp) + dfFacT*df(i,j)*maskUp(i,j)  
 c      ENDDO  
 c     ENDDO  
412    
413  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
414  C-    GM/Redi flux in R  C-    GM/Redi flux in R
# Line 304  C-    GM/Redi flux in R Line 416  C-    GM/Redi flux in R
416  C *note* should update GMREDI_RTRANSPORT to set df  *aja*  C *note* should update GMREDI_RTRANSPORT to set df  *aja*
417         CALL GMREDI_RTRANSPORT(         CALL GMREDI_RTRANSPORT(
418       I     iMin,iMax,jMin,jMax,bi,bj,K,       I     iMin,iMax,jMin,jMax,bi,bj,K,
419       I     maskUp,Tracer,       I     Tracer,tracerIdentity,
420       U     df,       U     df,
421       I     myThid)       I     myThid)
 c      DO j=jMin,jMax  
 c       DO i=iMin,iMax  
 c        fVerT(i,j,kUp) = fVerT(i,j,kUp) + dfFacT*df(i,j)*maskUp(i,j)  
 c       ENDDO  
 c      ENDDO  
422        ENDIF        ENDIF
423  #endif  #endif
424    
425        DO j=jMin,jMax        DO j=1-Oly,sNy+Oly
426         DO i=iMin,iMax         DO i=1-Olx,sNx+Olx
427          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)
428         ENDDO         ENDDO
429        ENDDO        ENDDO
430    
431    #ifdef ALLOW_DIAGNOSTICS
432    C-    Diagnostics of Tracer flux in R dir (mainly Diffusive terms),
433    C       Explicit terms only & excluding advective terms:
434          IF ( useDiagnostics .AND.
435         &    (.NOT.implicitDiffusion .OR. useGMRedi) ) THEN
436              diagName = 'DFrE'//diagSufx
437              CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
438          ENDIF
439    #endif
440    
441  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
442  C-    Add non local KPP transport term (ghat) to diffusive T flux.  C-    Set non local KPP transport term (ghat):
443        IF (useKPP) THEN        IF ( useKPP .AND. k.GE.2 ) THEN
444         DO j=jMin,jMax         DO j=1-Oly,sNy+Oly
445          DO i=iMin,iMax          DO i=1-Olx,sNx+Olx
446           df(i,j) = 0.           df(i,j) = 0. _d 0
447          ENDDO          ENDDO
448         ENDDO         ENDDO
449         IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN         IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
 C *note* should update KPP_TRANSPORT_T to set df  *aja*  
450          CALL KPP_TRANSPORT_T(          CALL KPP_TRANSPORT_T(
451       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,
452       I     KappaRT,       O     df )
      U     df )  
453         ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN         ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
454          CALL KPP_TRANSPORT_S(          CALL KPP_TRANSPORT_S(
455       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,
456       I     KappaRT,       O     df )
457       U     df )  #ifdef ALLOW_PTRACERS
458           ELSEIF (tracerIdentity .GE. GAD_TR1) THEN
459            CALL KPP_TRANSPORT_PTR(
460         I     iMin,iMax,jMin,jMax,bi,bj,k,km1,
461         I     tracerIdentity-GAD_TR1+1,
462         O     df )
463    #endif
464         ELSE         ELSE
465            PRINT*,'invalid tracer indentity: ', tracerIdentity
466          STOP 'GAD_CALC_RHS: Ooops'          STOP 'GAD_CALC_RHS: Ooops'
467         ENDIF         ENDIF
468         DO j=jMin,jMax         DO j=1-Oly,sNy+Oly
469          DO i=iMin,iMax          DO i=1-Olx,sNx+Olx
470           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)
471          ENDDO          ENDDO
472         ENDDO         ENDDO
473        ENDIF        ENDIF
474  #endif  #endif
475    
476  C--   Divergence of fluxes  C--   Divergence of fluxes
477        DO j=jMin,jMax        DO j=1-Oly,sNy+Oly-1
478         DO i=iMin,iMax         DO i=1-Olx,sNx+Olx-1
479          gTracer(i,j,k,bi,bj)=          gTracer(i,j,k,bi,bj)=gTracer(i,j,k,bi,bj)
480       &   -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &   -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)*recip_rA(i,j,bi,bj)
481       &    *recip_rA(i,j,bi,bj)       &   *( (fZon(i+1,j)-fZon(i,j))
482       &    *(       &     +(fMer(i,j+1)-fMer(i,j))
483       &    +( fZon(i+1,j)-fZon(i,j) )       &     +(fVerT(i,j,kUp)-fVerT(i,j,kDown))*rkFac
484       &    +( fMer(i,j+1)-fMer(i,j) )       &     -localT(i,j)*( (uTrans(i+1,j)-uTrans(i,j))
485       &    +( fVerT(i,j,kUp)-fVerT(i,j,kDown) )*rkFac       &                   +(vTrans(i,j+1)-vTrans(i,j))
486         &                   +(rTrans(i,j)-rTransKp1(i,j))*rAdvFac
487         &                  )*advFac
488       &    )       &    )
489         ENDDO         ENDDO
490        ENDDO        ENDDO
491    
492    #ifdef ALLOW_DEBUG
493          IF ( debugLevel .GE. debLevB
494         &   .AND. tracerIdentity.EQ.GAD_TEMPERATURE
495         &   .AND. k.EQ.2 .AND. myIter.EQ.1+nIter0
496         &   .AND. nPx.EQ.1 .AND. nPy.EQ.1
497         &   .AND. useCubedSphereExchange ) THEN
498            CALL DEBUG_CS_CORNER_UV( ' fZon,fMer from GAD_CALC_RHS',
499         &             fZon,fMer, k, standardMessageUnit,bi,bj,myThid )
500          ENDIF
501    #endif /* ALLOW_DEBUG */
502    
503        RETURN        RETURN
504        END        END
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.34
  ViewVC Help
Powered by ViewVC 1.1.22