/[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.12.6.1 by heimbach, Sun Mar 24 17:22:22 2002 UTC revision 1.38 by jmc, Sun Nov 6 22:14:02 2005 UTC
# Line 9  C !ROUTINE: GAD_CALC_RHS Line 9  C !ROUTINE: GAD_CALC_RHS
9  C !INTERFACE: ==========================================================  C !INTERFACE: ==========================================================
10        SUBROUTINE GAD_CALC_RHS(        SUBROUTINE GAD_CALC_RHS(
11       I           bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,       I           bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,
12       I           xA,yA,uTrans,vTrans,rTrans,maskUp,       I           xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,
13       I           diffKh, diffK4, KappaRT, Tracer,       I           uVel, vVel, wVel,
14       I           tracerIdentity, advectionScheme,       I           diffKh, diffK4, KappaR, Tracer, TracAB,
15         I           tracerIdentity, advectionScheme, vertAdvecScheme,
16         I           calcAdvection, implicitAdvection,
17       U           fVerT, gTracer,       U           fVerT, gTracer,
18       I           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 27  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 40  C !USES: =============================== Line 42  C !USES: ===============================
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  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 49  C !USES: =============================== Line 51  C !USES: ===============================
51  #endif /* ALLOW_AUTODIFF_TAMC */  #endif /* ALLOW_AUTODIFF_TAMC */
52    
53  C !INPUT PARAMETERS: ===================================================  C !INPUT PARAMETERS: ===================================================
54  C  bi,bj                :: tile indices  C bi,bj            :: tile indices
55  C  iMin,iMax,jMin,jMax  :: loop range for called routines  C iMin,iMax        :: loop range for called routines
56  C  kup                  :: index into 2 1/2D array, toggles between 1 and 2  C jMin,jMax        :: loop range for called routines
57  C  kdown                :: index into 2 1/2D array, toggles between 2 and 1  C kup              :: index into 2 1/2D array, toggles between 1|2
58  C  kp1                  :: =k+1 for k<Nr, =Nr for k=Nr  C kdown            :: index into 2 1/2D array, toggles between 2|1
59  C  xA,yA                :: areas of X and Y face of tracer cells  C kp1              :: =k+1 for k<Nr, =Nr for k=Nr
60  C  uTrans,vTrans,rTrans :: 2-D arrays of volume transports at U,V and W points  C xA,yA            :: areas of X and Y face of tracer cells
61  C  maskUp               :: 2-D array for mask at W points  C uTrans,vTrans    :: 2-D arrays of volume transports at U,V points
62  C  diffKh               :: horizontal diffusion coefficient  C rTrans           :: 2-D arrays of volume transports at W points
63  C  diffK4               :: bi-harmonic diffusion coefficient  C rTransKp1        :: 2-D array of volume trans at W pts, interf k+1
64  C  KappaRT              :: 3-D array for vertical diffusion coefficient  C maskUp           :: 2-D array for mask at W points
65  C  Tracer               :: tracer field  C uVel,vVel,wVel   :: 3 components of the velcity field (3-D array)
66  C  tracerIdentity       :: identifier for the tracer (required only for KPP)  C diffKh           :: horizontal diffusion coefficient
67  C  advectionScheme      :: advection scheme to use  C diffK4           :: bi-harmonic diffusion coefficient
68  C  myThid               :: thread number  C KappaR           :: 2-D array for vertical diffusion coefficient, interf k
69    C Tracer           :: tracer field
70    C TracAB           :: tracer field but extrapolated fwd in time (to nIter+1/2)
71    C                     if applying AB on tracer field (instead of on tendency)
72    C tracerIdentity   :: tracer identifier (required for KPP,GM)
73    C advectionScheme  :: advection scheme to use (Horizontal plane)
74    C vertAdvecScheme  :: advection scheme to use (Vertical direction)
75    C calcAdvection    :: =False if Advec computed with multiDim scheme
76    C implicitAdvection:: =True if vertical Advec computed implicitly
77    C myTime           :: current time
78    C myIter           :: iteration number
79    C myThid           :: thread number
80        INTEGER bi,bj,iMin,iMax,jMin,jMax        INTEGER bi,bj,iMin,iMax,jMin,jMax
81        INTEGER k,kUp,kDown,kM1        INTEGER k,kUp,kDown,kM1
82        _RS xA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS xA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 71  C  myThid               :: thread number Line 84  C  myThid               :: thread number
84        _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
85        _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
86        _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
87          _RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
88        _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
89          _RL uVel  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
90          _RL vVel  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
91          _RL wVel  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
92        _RL diffKh, diffK4        _RL diffKh, diffK4
93        _RL KappaRT(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL KappaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
94        _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)
95          _RL TracAB(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
96        INTEGER tracerIdentity        INTEGER tracerIdentity
97        INTEGER advectionScheme        INTEGER advectionScheme, vertAdvecScheme
98        INTEGER myThid        LOGICAL calcAdvection
99          LOGICAL implicitAdvection
100          _RL     myTime
101          INTEGER myIter, myThid
102    
103  C !OUTPUT PARAMETERS: ==================================================  C !OUTPUT PARAMETERS: ==================================================
104  C  gTracer              :: tendancy array  C gTracer          :: tendency array
105  C  fVerT                :: 2 1/2D arrays for vertical advective flux  C fVerT            :: 2 1/2D arrays for vertical advective flux
106        _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)
107        _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
108    
109  C !LOCAL VARIABLES: ====================================================  C !LOCAL VARIABLES: ====================================================
110  C  i,j                  :: loop indices  C i,j              :: loop indices
111  C  df4                  :: used for storing del^2 T for bi-harmonic term  C df4              :: used for storing del^2 T for bi-harmonic term
112  C  fZon                 :: zonal flux  C fZon             :: zonal flux
113  C  fmer                 :: meridional flux  C fMer             :: meridional flux
114  C  af                   :: advective flux  C af               :: advective flux
115  C  df                   :: diffusive flux  C df               :: diffusive flux
116  C  localT               :: local copy of tracer field  C localT           :: local copy of tracer field
117    C locABT           :: local copy of (AB-extrapolated) tracer field
118    #ifdef ALLOW_DIAGNOSTICS
119          CHARACTER*8 diagName
120          CHARACTER*4 GAD_DIAG_SUFX, diagSufx
121          EXTERNAL    GAD_DIAG_SUFX
122    #endif
123        INTEGER i,j        INTEGER i,j
124        _RL df4   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL df4   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
125        _RL fZon  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fZon  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 100  C  localT               :: local copy of Line 127  C  localT               :: local copy of
127        _RL af    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL af    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
128        _RL df    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL df    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
129        _RL localT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL localT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
130          _RL locABT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
131          _RL advFac, rAdvFac
132  CEOP  CEOP
133    
134  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 108  C--   the kDown is still required Line 137  C--   the kDown is still required
137        fVerT(1,1,kDown) = fVerT(1,1,kDown)        fVerT(1,1,kDown) = fVerT(1,1,kDown)
138  #endif  #endif
139    
140    #ifdef ALLOW_DIAGNOSTICS
141    C--   Set diagnostic suffix for the current tracer
142          IF ( useDiagnostics ) THEN
143            diagSufx = GAD_DIAG_SUFX( tracerIdentity, myThid )
144          ENDIF
145    #endif
146    
147          advFac  = 0. _d 0
148          IF (calcAdvection) advFac = 1. _d 0
149          rAdvFac = rkSign*advFac
150          IF (implicitAdvection) rAdvFac = 0. _d 0
151    
152        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
153         DO i=1-OLx,sNx+OLx         DO i=1-OLx,sNx+OLx
154          fZon(i,j)      = 0. _d 0          fZon(i,j)      = 0. _d 0
# Line 115  C--   the kDown is still required Line 156  C--   the kDown is still required
156          fVerT(i,j,kUp) = 0. _d 0          fVerT(i,j,kUp) = 0. _d 0
157          df(i,j)        = 0. _d 0          df(i,j)        = 0. _d 0
158          df4(i,j)       = 0. _d 0          df4(i,j)       = 0. _d 0
         localT(i,j)    = 0. _d 0  
159         ENDDO         ENDDO
160        ENDDO        ENDDO
161    
162  C--   Make local copy of tracer array  C--   Make local copy of tracer array
163        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
164         DO i=1-OLx,sNx+OLx         DO i=1-OLx,sNx+OLx
165          localT(i,j)=tracer(i,j,k,bi,bj)          localT(i,j)=Tracer(i,j,k,bi,bj)
166            locABT(i,j)=TracAB(i,j,k,bi,bj)
167         ENDDO         ENDDO
168        ENDDO        ENDDO
169    
170  C--   Unless we have already calculated the advection terms we initialize  C--   Unless we have already calculated the advection terms we initialize
171  C     the tendency to zero.  C     the tendency to zero.
172        IF (.NOT. multiDimAdvection .OR.  C     <== now done earlier at the beginning of thermodynamics.
173       &    advectionScheme.EQ.ENUM_CENTERED_2ND .OR.  c     IF (calcAdvection) THEN
174       &    advectionScheme.EQ.ENUM_UPWIND_3RD .OR.  c      DO j=1-Oly,sNy+Oly
175       &    advectionScheme.EQ.ENUM_CENTERED_4TH ) THEN  c       DO i=1-Olx,sNx+Olx
176         DO j=1-Oly,sNy+Oly  c        gTracer(i,j,k,bi,bj)=0. _d 0
177          DO i=1-Olx,sNx+Olx  c       ENDDO
178           gTracer(i,j,k,bi,bj)=0. _d 0  c      ENDDO
179          ENDDO  c     ENDIF
        ENDDO  
       ENDIF  
180    
181  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
182        IF (diffK4 .NE. 0.) THEN        IF (diffK4 .NE. 0.) THEN
# Line 154  C--   Initialize net flux in X direction Line 193  C--   Initialize net flux in X direction
193        ENDDO        ENDDO
194    
195  C-    Advective flux in X  C-    Advective flux in X
196        IF (.NOT. multiDimAdvection .OR.        IF (calcAdvection) THEN
197       &    advectionScheme.EQ.ENUM_CENTERED_2ND .OR.          IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
198       &    advectionScheme.EQ.ENUM_UPWIND_3RD .OR.            CALL GAD_C2_ADV_X(bi,bj,k,uTrans,locABT,af,myThid)
199       &    advectionScheme.EQ.ENUM_CENTERED_4TH ) THEN          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST
200        IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN       &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
201         CALL GAD_C2_ADV_X(bi,bj,k,uTrans,localT,af,myThid)            CALL GAD_DST2U1_ADV_X( bi,bj,k, advectionScheme,
202        ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN       I            dTtracerLev(k), uTrans, uVel, locABT,
203         CALL GAD_FLUXLIMIT_ADV_X(       O            af, myThid )
204       &      bi,bj,k,deltaTtracer,uTrans,uVel,localT,af,myThid)          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
205        ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN            CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, dTtracerLev(k),
206         CALL GAD_U3_ADV_X(bi,bj,k,uTrans,localT,af,myThid)       I            uTrans, uVel, maskW(1-Olx,1-Oly,k,bi,bj), locABT,
207        ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN       O            af, myThid )
208         CALL GAD_C4_ADV_X(bi,bj,k,uTrans,localT,af,myThid)          ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
209        ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN            CALL GAD_U3_ADV_X(bi,bj,k,uTrans,locABT,af,myThid)
210         CALL GAD_DST3_ADV_X(          ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
211       &       bi,bj,k,deltaTtracer,uTrans,uVel,localT,af,myThid)            CALL GAD_C4_ADV_X(bi,bj,k,uTrans,locABT,af,myThid)
212        ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
213         CALL GAD_DST3FL_ADV_X(            CALL GAD_DST3_ADV_X( bi,bj,k, dTtracerLev(k),
214       &       bi,bj,k,deltaTtracer,uTrans,uVel,localT,af,myThid)       I            uTrans, uVel, maskW(1-Olx,1-Oly,k,bi,bj), locABT,
215        ELSE       O            af, myThid )
216         STOP 'GAD_CALC_RHS: Bad advectionScheme (X)'          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
217        ENDIF           IF ( inAdMode ) THEN
218        DO j=1-Oly,sNy+Oly  cph This block is to trick the adjoint:
219         DO i=1-Olx,sNx+Olx  cph IF inAdExact=.FALSE., we want to use DST3
220          fZon(i,j) = fZon(i,j) + af(i,j)  cph with limiters in forward, but without limiters in reverse.
221         ENDDO            CALL GAD_DST3_ADV_X( bi,bj,k, dTtracerLev(k),
222        ENDDO       I           uTrans, uVel, maskW(1-Olx,1-Oly,k,bi,bj), locABT,
223         O           af, myThid )
224             ELSE
225              CALL GAD_DST3FL_ADV_X( bi,bj,k, dTtracerLev(k),
226         I           uTrans, uVel, maskW(1-Olx,1-Oly,k,bi,bj), locABT,
227         O           af, myThid )
228             ENDIF
229            ELSE
230             STOP 'GAD_CALC_RHS: Bad advectionScheme (X)'
231            ENDIF
232            DO j=1-Oly,sNy+Oly
233             DO i=1-Olx,sNx+Olx
234              fZon(i,j) = fZon(i,j) + af(i,j)
235             ENDDO
236            ENDDO
237    #ifdef ALLOW_DIAGNOSTICS
238            IF ( useDiagnostics ) THEN
239              diagName = 'ADVx'//diagSufx
240              CALL DIAGNOSTICS_FILL(af,diagName, k,1, 2,bi,bj, myThid)
241            ENDIF
242    #endif
243        ENDIF        ENDIF
244    
245  C-    Diffusive flux in X  C-    Diffusive flux in X
# Line 194  C-    Diffusive flux in X Line 253  C-    Diffusive flux in X
253         ENDDO         ENDDO
254        ENDIF        ENDIF
255    
256    C-    Add bi-harmonic diffusive flux in X
257          IF (diffK4 .NE. 0.) THEN
258           CALL GAD_BIHARM_X(bi,bj,k,xA,df4,diffK4,df,myThid)
259          ENDIF
260    
261  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
262  C-    GM/Redi flux in X  C-    GM/Redi flux in X
263        IF (useGMRedi) THEN        IF (useGMRedi) THEN
264  C *note* should update GMREDI_XTRANSPORT to use localT and set df  *aja*  C *note* should update GMREDI_XTRANSPORT to set df  *aja*
265          CALL GMREDI_XTRANSPORT(          CALL GMREDI_XTRANSPORT(
266       I     iMin,iMax,jMin,jMax,bi,bj,K,       I     iMin,iMax,jMin,jMax,bi,bj,K,
267       I     xA,Tracer,       I     xA,Tracer,tracerIdentity,
268       U     df,       U     df,
269       I     myThid)       I     myThid)
270        ENDIF        ENDIF
# Line 211  C *note* should update GMREDI_XTRANSPORT Line 275  C *note* should update GMREDI_XTRANSPORT
275         ENDDO         ENDDO
276        ENDDO        ENDDO
277    
278  C-    Bi-harmonic duffusive flux in X  #ifdef ALLOW_DIAGNOSTICS
279        IF (diffK4 .NE. 0.) THEN  C-    Diagnostics of Tracer flux in X dir (mainly Diffusive term),
280         CALL GAD_BIHARM_X(bi,bj,k,xA,df4,diffK4,df,myThid)  C       excluding advective terms:
281         DO j=1-Oly,sNy+Oly        IF ( useDiagnostics .AND.
282          DO i=1-Olx,sNx+Olx       &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. useGMRedi) ) THEN
283           fZon(i,j) = fZon(i,j) + df(i,j)            diagName = 'DIFx'//diagSufx
284          ENDDO            CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
        ENDDO  
285        ENDIF        ENDIF
286    #endif
287    
288  C--   Initialize net flux in Y direction  C--   Initialize net flux in Y direction
289        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
# Line 229  C--   Initialize net flux in Y direction Line 293  C--   Initialize net flux in Y direction
293        ENDDO        ENDDO
294    
295  C-    Advective flux in Y  C-    Advective flux in Y
296        IF (.NOT. multiDimAdvection .OR.        IF (calcAdvection) THEN
297       &    advectionScheme.EQ.ENUM_CENTERED_2ND .OR.          IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
298       &    advectionScheme.EQ.ENUM_UPWIND_3RD .OR.            CALL GAD_C2_ADV_Y(bi,bj,k,vTrans,locABT,af,myThid)
299       &    advectionScheme.EQ.ENUM_CENTERED_4TH ) THEN          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST
300        IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN       &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
301         CALL GAD_C2_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)            CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme,
302        ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN       I            dTtracerLev(k), vTrans, vVel, locABT,
303         CALL GAD_FLUXLIMIT_ADV_Y(       O            af, myThid )
304       &       bi,bj,k,deltaTtracer,vTrans,vVel,localT,af,myThid)          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
305        ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN            CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, dTtracerLev(k),
306         CALL GAD_U3_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)       I            vTrans, vVel, maskS(1-Olx,1-Oly,k,bi,bj), locABT,
307        ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN       O            af, myThid )
308         CALL GAD_C4_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)          ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
309        ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN            CALL GAD_U3_ADV_Y(bi,bj,k,vTrans,locABT,af,myThid)
310         CALL GAD_DST3_ADV_Y(          ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
311       &       bi,bj,k,deltaTtracer,vTrans,vVel,localT,af,myThid)            CALL GAD_C4_ADV_Y(bi,bj,k,vTrans,locABT,af,myThid)
312        ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
313         CALL GAD_DST3FL_ADV_Y(            CALL GAD_DST3_ADV_Y( bi,bj,k, dTtracerLev(k),
314       &       bi,bj,k,deltaTtracer,vTrans,vVel,localT,af,myThid)       I            vTrans, vVel, maskS(1-Olx,1-Oly,k,bi,bj), locABT,
315        ELSE       O            af, myThid )
316         STOP 'GAD_CALC_RHS: Bad advectionScheme (Y)'          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
317        ENDIF           IF ( inAdMode ) THEN
318        DO j=1-Oly,sNy+Oly  cph This block is to trick the adjoint:
319         DO i=1-Olx,sNx+Olx  cph IF inAdExact=.FALSE., we want to use DST3
320          fMer(i,j) = fMer(i,j) + af(i,j)  cph with limiters in forward, but without limiters in reverse.
321         ENDDO            CALL GAD_DST3_ADV_Y( bi,bj,k, dTtracerLev(k),
322        ENDDO       I           vTrans, vVel, maskS(1-Olx,1-Oly,k,bi,bj), locABT,
323         O           af, myThid )
324             ELSE
325              CALL GAD_DST3FL_ADV_Y( bi,bj,k, dTtracerLev(k),
326         I           vTrans, vVel, maskS(1-Olx,1-Oly,k,bi,bj), locABT,
327         O           af, myThid )
328             ENDIF
329            ELSE
330              STOP 'GAD_CALC_RHS: Bad advectionScheme (Y)'
331            ENDIF
332            DO j=1-Oly,sNy+Oly
333             DO i=1-Olx,sNx+Olx
334              fMer(i,j) = fMer(i,j) + af(i,j)
335             ENDDO
336            ENDDO
337    #ifdef ALLOW_DIAGNOSTICS
338            IF ( useDiagnostics ) THEN
339              diagName = 'ADVy'//diagSufx
340              CALL DIAGNOSTICS_FILL(af,diagName, k,1, 2,bi,bj, myThid)
341            ENDIF
342    #endif
343        ENDIF        ENDIF
344    
345  C-    Diffusive flux in Y  C-    Diffusive flux in Y
# Line 269  C-    Diffusive flux in Y Line 353  C-    Diffusive flux in Y
353         ENDDO         ENDDO
354        ENDIF        ENDIF
355    
356    C-    Add bi-harmonic flux in Y
357          IF (diffK4 .NE. 0.) THEN
358           CALL GAD_BIHARM_Y(bi,bj,k,yA,df4,diffK4,df,myThid)
359          ENDIF
360    
361  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
362  C-    GM/Redi flux in Y  C-    GM/Redi flux in Y
363        IF (useGMRedi) THEN        IF (useGMRedi) THEN
364  C *note* should update GMREDI_YTRANSPORT to use localT and set df  *aja*  C *note* should update GMREDI_YTRANSPORT to set df  *aja*
365         CALL GMREDI_YTRANSPORT(         CALL GMREDI_YTRANSPORT(
366       I     iMin,iMax,jMin,jMax,bi,bj,K,       I     iMin,iMax,jMin,jMax,bi,bj,K,
367       I     yA,Tracer,       I     yA,Tracer,tracerIdentity,
368       U     df,       U     df,
369       I     myThid)       I     myThid)
370        ENDIF        ENDIF
# Line 286  C *note* should update GMREDI_YTRANSPORT Line 375  C *note* should update GMREDI_YTRANSPORT
375         ENDDO         ENDDO
376        ENDDO        ENDDO
377    
378  C-    Bi-harmonic flux in Y  #ifdef ALLOW_DIAGNOSTICS
379        IF (diffK4 .NE. 0.) THEN  C-    Diagnostics of Tracer flux in Y dir (mainly Diffusive terms),
380         CALL GAD_BIHARM_Y(bi,bj,k,yA,df4,diffK4,df,myThid)  C       excluding advective terms:
381         DO j=1-Oly,sNy+Oly        IF ( useDiagnostics .AND.
382          DO i=1-Olx,sNx+Olx       &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. useGMRedi) ) THEN
383           fMer(i,j) = fMer(i,j) + df(i,j)            diagName = 'DIFy'//diagSufx
384          ENDDO            CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
        ENDDO  
385        ENDIF        ENDIF
386    #endif
387    
388    C--   Compute vertical flux fVerT(kUp) at interface k (between k-1 & k):
389  C-    Advective flux in R  C-    Advective flux in R
390        IF (.NOT. multiDimAdvection .OR.  #ifdef ALLOW_AIM
391       &    advectionScheme.EQ.ENUM_CENTERED_2ND .OR.  C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr
392       &    advectionScheme.EQ.ENUM_UPWIND_3RD .OR.        IF (calcAdvection .AND. .NOT.implicitAdvection .AND. K.GE.2 .AND.
393       &    advectionScheme.EQ.ENUM_CENTERED_4TH ) THEN       &     (.NOT.useAIM .OR.tracerIdentity.NE.GAD_SALINITY .OR.K.LT.Nr)
394  C     Note: wVel needs to be masked       &   ) THEN
395        IF (K.GE.2) THEN  #else
396          IF (calcAdvection .AND. .NOT.implicitAdvection .AND. K.GE.2) THEN
397    #endif
398  C-    Compute vertical advective flux in the interior:  C-    Compute vertical advective flux in the interior:
399         IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN          IF (vertAdvecScheme.EQ.ENUM_CENTERED_2ND) THEN
400          CALL GAD_C2_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)            CALL GAD_C2_ADV_R(bi,bj,k,rTrans,TracAB,af,myThid)
401         ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF ( vertAdvecScheme.EQ.ENUM_UPWIND_1RST
402          CALL GAD_FLUXLIMIT_ADV_R(       &          .OR. vertAdvecScheme.EQ.ENUM_DST2 ) THEN
403       &       bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)            CALL GAD_DST2U1_ADV_R( bi,bj,k, vertAdvecScheme,
404         ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN       I         dTtracerLev(k),rTrans,wVel,TracAB(1-Olx,1-Oly,1,bi,bj),
405          CALL GAD_U3_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)       O         af, myThid )
406         ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN          ELSEIF (vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN
407          CALL GAD_C4_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)            CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k,
408         ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN       I         dTtracerLev(k),rTrans,wVel,TracAB(1-Olx,1-Oly,1,bi,bj),
409          CALL GAD_DST3_ADV_R(       O         af, myThid )
410       &       bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)          ELSEIF (vertAdvecScheme.EQ.ENUM_UPWIND_3RD ) THEN
411         ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN            CALL GAD_U3_ADV_R(bi,bj,k,rTrans,TracAB,af,myThid)
412          CALL GAD_DST3FL_ADV_R(          ELSEIF (vertAdvecScheme.EQ.ENUM_CENTERED_4TH) THEN
413       &       bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)            CALL GAD_C4_ADV_R(bi,bj,k,rTrans,TracAB,af,myThid)
414         ELSE          ELSEIF (vertAdvecScheme.EQ.ENUM_DST3 ) THEN
415          STOP 'GAD_CALC_RHS: Bad advectionScheme (R)'            CALL GAD_DST3_ADV_R( bi,bj,k,
416         ENDIF       I         dTtracerLev(k),rTrans,wVel,TracAB(1-Olx,1-Oly,1,bi,bj),
417  C-    Surface "correction" term at k>1 :       O         af, myThid )
418         DO j=1-Oly,sNy+Oly          ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
419          DO i=1-Olx,sNx+Olx  cph This block is to trick the adjoint:
420           af(i,j) = af(i,j)  cph IF inAdExact=.FALSE., we want to use DST3
421       &           + (maskC(i,j,k,bi,bj)-maskC(i,j,k-1,bi,bj))*  cph with limiters in forward, but without limiters in reverse.
422       &             rTrans(i,j)*Tracer(i,j,k,bi,bj)            IF ( inAdMode ) THEN
423               CALL GAD_DST3_ADV_R( bi,bj,k,
424         I         dTtracerLev(k),rTrans,wVel,TracAB(1-Olx,1-Oly,1,bi,bj),
425         O         af, myThid )
426              ELSE
427               CALL GAD_DST3FL_ADV_R( bi,bj,k,
428         I         dTtracerLev(k),rTrans,wVel,TracAB(1-Olx,1-Oly,1,bi,bj),
429         O         af, myThid )
430              ENDIF
431            ELSE
432              STOP 'GAD_CALC_RHS: Bad vertAdvecScheme (R)'
433            ENDIF
434    C-     add the advective flux to fVerT
435            DO j=1-Oly,sNy+Oly
436             DO i=1-Olx,sNx+Olx
437              fVerT(i,j,kUp) = fVerT(i,j,kUp) + af(i,j)
438             ENDDO
439          ENDDO          ENDDO
440         ENDDO  #ifdef ALLOW_DIAGNOSTICS
441        ELSE          IF ( useDiagnostics ) THEN
442  C-    Surface "correction" term at k=1 :            diagName = 'ADVr'//diagSufx
443         DO j=1-Oly,sNy+Oly            CALL DIAGNOSTICS_FILL(af,diagName, k,1, 2,bi,bj, myThid)
444          DO i=1-Olx,sNx+Olx  C- note: needs to explicitly increment the counter since DIAGNOSTICS_FILL
445           af(i,j) = rTrans(i,j)*Tracer(i,j,k,bi,bj)  C        does it only if k=1 (never the case here)
446          ENDDO            IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT(diagName,bi,bj,myThid)
447         ENDDO          ENDIF
448        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) + af(i,j)  
        ENDDO  
       ENDDO  
449        ENDIF        ENDIF
450    
451  C-    Diffusive flux in R  C-    Diffusive flux in R
# Line 356  C           boundary condition. Line 458  C           boundary condition.
458          ENDDO          ENDDO
459         ENDDO         ENDDO
460        ELSE        ELSE
461         CALL GAD_DIFF_R(bi,bj,k,KappaRT,tracer,df,myThid)         CALL GAD_DIFF_R(bi,bj,k,KappaR,Tracer,df,myThid)
462        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) + df(i,j)*maskUp(i,j)  
 c      ENDDO  
 c     ENDDO  
463    
464  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
465  C-    GM/Redi flux in R  C-    GM/Redi flux in R
# Line 370  C-    GM/Redi flux in R Line 467  C-    GM/Redi flux in R
467  C *note* should update GMREDI_RTRANSPORT to set df  *aja*  C *note* should update GMREDI_RTRANSPORT to set df  *aja*
468         CALL GMREDI_RTRANSPORT(         CALL GMREDI_RTRANSPORT(
469       I     iMin,iMax,jMin,jMax,bi,bj,K,       I     iMin,iMax,jMin,jMax,bi,bj,K,
470       I     Tracer,       I     Tracer,tracerIdentity,
471       U     df,       U     df,
472       I     myThid)       I     myThid)
 c      DO j=1-Oly,sNy+Oly  
 c       DO i=1-Olx,sNx+Olx  
 c        fVerT(i,j,kUp) = fVerT(i,j,kUp) + df(i,j)*maskUp(i,j)  
 c       ENDDO  
 c      ENDDO  
473        ENDIF        ENDIF
474  #endif  #endif
475    
# Line 387  c      ENDDO Line 479  c      ENDDO
479         ENDDO         ENDDO
480        ENDDO        ENDDO
481    
482    #ifdef ALLOW_DIAGNOSTICS
483    C-    Diagnostics of Tracer flux in R dir (mainly Diffusive terms),
484    C       Explicit terms only & excluding advective terms:
485          IF ( useDiagnostics .AND.
486         &    (.NOT.implicitDiffusion .OR. useGMRedi) ) THEN
487              diagName = 'DFrE'//diagSufx
488              CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
489          ENDIF
490    #endif
491    
492  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
493  C-    Add non local KPP transport term (ghat) to diffusive T flux.  C-    Set non local KPP transport term (ghat):
494        IF (useKPP) THEN        IF ( useKPP .AND. k.GE.2 ) THEN
495         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
496          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
497           df(i,j) = 0. _d 0           df(i,j) = 0. _d 0
498          ENDDO          ENDDO
499         ENDDO         ENDDO
500         IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN         IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
 C *note* should update KPP_TRANSPORT_T to set df  *aja*  
501          CALL KPP_TRANSPORT_T(          CALL KPP_TRANSPORT_T(
502       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,
503       I     KappaRT,       O     df )
      U     df )  
504         ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN         ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
505          CALL KPP_TRANSPORT_S(          CALL KPP_TRANSPORT_S(
506       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,
507       I     KappaRT,       O     df )
508       U     df )  #ifdef ALLOW_PTRACERS
509           ELSEIF (tracerIdentity .GE. GAD_TR1) THEN
510            CALL KPP_TRANSPORT_PTR(
511         I     iMin,iMax,jMin,jMax,bi,bj,k,km1,
512         I     tracerIdentity-GAD_TR1+1,
513         O     df )
514    #endif
515         ELSE         ELSE
516            PRINT*,'invalid tracer indentity: ', tracerIdentity
517          STOP 'GAD_CALC_RHS: Ooops'          STOP 'GAD_CALC_RHS: Ooops'
518         ENDIF         ENDIF
519         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
# Line 421  C--   Divergence of fluxes Line 528  C--   Divergence of fluxes
528        DO j=1-Oly,sNy+Oly-1        DO j=1-Oly,sNy+Oly-1
529         DO i=1-Olx,sNx+Olx-1         DO i=1-Olx,sNx+Olx-1
530          gTracer(i,j,k,bi,bj)=gTracer(i,j,k,bi,bj)          gTracer(i,j,k,bi,bj)=gTracer(i,j,k,bi,bj)
531       &   -_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)
532       &    *recip_rA(i,j,bi,bj)       &   *( (fZon(i+1,j)-fZon(i,j))
533       &    *(       &     +(fMer(i,j+1)-fMer(i,j))
534       &    +( fZon(i+1,j)-fZon(i,j) )       &     +(fVerT(i,j,kDown)-fVerT(i,j,kUp))*rkSign
535       &    +( fMer(i,j+1)-fMer(i,j) )       &     -localT(i,j)*( (uTrans(i+1,j)-uTrans(i,j))
536       &    +( fVerT(i,j,kUp)-fVerT(i,j,kDown) )*rkFac       &                   +(vTrans(i,j+1)-vTrans(i,j))
537         &                   +(rTransKp1(i,j)-rTrans(i,j))*rAdvFac
538         &                  )*advFac
539       &    )       &    )
540         ENDDO         ENDDO
541        ENDDO        ENDDO
542    
543    #ifdef ALLOW_DEBUG
544          IF ( debugLevel .GE. debLevB
545         &   .AND. tracerIdentity.EQ.GAD_TEMPERATURE
546         &   .AND. k.EQ.2 .AND. myIter.EQ.1+nIter0
547         &   .AND. nPx.EQ.1 .AND. nPy.EQ.1
548         &   .AND. useCubedSphereExchange ) THEN
549            CALL DEBUG_CS_CORNER_UV( ' fZon,fMer from GAD_CALC_RHS',
550         &             fZon,fMer, k, standardMessageUnit,bi,bj,myThid )
551          ENDIF
552    #endif /* ALLOW_DEBUG */
553    
554        RETURN        RETURN
555        END        END
Legend:
Removed from v.1.12.6.1  
changed lines
  Added in v.1.38
  ViewVC Help
Powered by ViewVC 1.1.22