/[MITgcm]/MITgcm/pkg/generic_advdiff/gad_calc_rhs.F
ViewVC logotype

Diff of /MITgcm/pkg/generic_advdiff/gad_calc_rhs.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.17 by jmc, Sun Jan 26 21:08:36 2003 UTC revision 1.37 by jmc, Sat Oct 22 20:09:52 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, calcAdvection,       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    
20  C !DESCRIPTION:  C !DESCRIPTION:
21  C Calculates the tendancy of a tracer due to advection and diffusion.  C Calculates the tendancy of a tracer due to advection and diffusion.
# 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"
 #include "DYNVARS.h"  
45  #include "SURFACE.h"  #include "SURFACE.h"
46  #include "GAD.h"  #include "GAD.h"
47    
# Line 50  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  calcAdvection        :: =False if Advec terms computed with multiDim scheme  C KappaR           :: 2-D array for vertical diffusion coefficient, interf k
69  C  myThid               :: thread number  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        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)
# Line 73  C  myThid               :: thread number Line 82  C  myThid               :: thread number
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        LOGICAL calcAdvection        LOGICAL calcAdvection
96        INTEGER myThid        LOGICAL implicitAdvection
97          _RL     myTime
98          INTEGER myIter, myThid
99    
100  C !OUTPUT PARAMETERS: ==================================================  C !OUTPUT PARAMETERS: ==================================================
101  C  gTracer              :: tendancy array  C gTracer          :: tendancy array
102  C  fVerT                :: 2 1/2D arrays for vertical advective flux  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        _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
105    
106  C !LOCAL VARIABLES: ====================================================  C !LOCAL VARIABLES: ====================================================
107  C  i,j                  :: loop indices  C i,j              :: loop indices
108  C  df4                  :: used for storing del^2 T for bi-harmonic term  C df4              :: used for storing del^2 T for bi-harmonic term
109  C  fZon                 :: zonal flux  C fZon             :: zonal flux
110  C  fmer                 :: meridional flux  C fMer             :: meridional flux
111  C  af                   :: advective flux  C af               :: advective flux
112  C  df                   :: diffusive flux  C df               :: diffusive flux
113  C  localT               :: local copy of tracer field  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
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)
# Line 103  C  localT               :: local copy of Line 123  C  localT               :: local copy of
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  CEOP
128    
129  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 111  C--   the kDown is still required Line 132  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 = rkSign*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. _d 0          fZon(i,j)      = 0. _d 0
# Line 118  C--   the kDown is still required Line 151  C--   the kDown is still required
151          fVerT(i,j,kUp) = 0. _d 0          fVerT(i,j,kUp) = 0. _d 0
152          df(i,j)        = 0. _d 0          df(i,j)        = 0. _d 0
153          df4(i,j)       = 0. _d 0          df4(i,j)       = 0. _d 0
         localT(i,j)    = 0. _d 0  
154         ENDDO         ENDDO
155        ENDDO        ENDDO
156    
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
160          localT(i,j)=tracer(i,j,k,bi,bj)          localT(i,j)=Tracer(i,j,k,bi,bj)
161         ENDDO         ENDDO
162        ENDDO        ENDDO
163    
164  C--   Unless we have already calculated the advection terms we initialize  C--   Unless we have already calculated the advection terms we initialize
165  C     the tendency to zero.  C     the tendency to zero.
166        IF (calcAdvection) THEN  C     <== now done earlier at the beginning of thermodynamics.
167         DO j=1-Oly,sNy+Oly  c     IF (calcAdvection) THEN
168          DO i=1-Olx,sNx+Olx  c      DO j=1-Oly,sNy+Oly
169           gTracer(i,j,k,bi,bj)=0. _d 0  c       DO i=1-Olx,sNx+Olx
170          ENDDO  c        gTracer(i,j,k,bi,bj)=0. _d 0
171         ENDDO  c       ENDDO
172        ENDIF  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 155  C--   Initialize net flux in X direction Line 188  C--   Initialize net flux in X direction
188    
189  C-    Advective flux in X  C-    Advective flux in X
190        IF (calcAdvection) THEN        IF (calcAdvection) THEN
191        IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN          IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
192         CALL GAD_C2_ADV_X(bi,bj,k,uTrans,localT,af,myThid)            CALL GAD_C2_ADV_X(bi,bj,k,uTrans,localT,af,myThid)
193        ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST
194         CALL GAD_FLUXLIMIT_ADV_X(       &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
195       &      bi,bj,k,deltaTtracer,uTrans,uVel,localT,af,myThid)            CALL GAD_DST2U1_ADV_X( bi,bj,k, advectionScheme,
196        ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN       I            dTtracerLev(k), uTrans, uVel, localT,
197         CALL GAD_U3_ADV_X(bi,bj,k,uTrans,localT,af,myThid)       O            af, myThid )
198        ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
199         CALL GAD_C4_ADV_X(bi,bj,k,uTrans,localT,af,myThid)            CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, dTtracerLev(k),
200        ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN       I            uTrans, uVel, maskW(1-Olx,1-Oly,k,bi,bj), localT,
201         CALL GAD_DST3_ADV_X(       O            af, myThid )
202       &       bi,bj,k,deltaTtracer,uTrans,uVel,localT,af,myThid)          ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
203        ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN            CALL GAD_U3_ADV_X(bi,bj,k,uTrans,localT,af,myThid)
204         CALL GAD_DST3FL_ADV_X(          ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
205       &       bi,bj,k,deltaTtracer,uTrans,uVel,localT,af,myThid)            CALL GAD_C4_ADV_X(bi,bj,k,uTrans,localT,af,myThid)
206        ELSE          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
207         STOP 'GAD_CALC_RHS: Bad advectionScheme (X)'            CALL GAD_DST3_ADV_X( bi,bj,k, dTtracerLev(k),
208        ENDIF       I            uTrans, uVel, maskW(1-Olx,1-Oly,k,bi,bj), localT,
209        DO j=1-Oly,sNy+Oly       O            af, myThid )
210         DO i=1-Olx,sNx+Olx          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
211          fZon(i,j) = fZon(i,j) + af(i,j)           IF ( inAdMode ) THEN
212         ENDDO  cph This block is to trick the adjoint:
213        ENDDO  cph IF inAdExact=.FALSE., we want to use DST3
214    cph with limiters in forward, but without limiters in reverse.
215              CALL GAD_DST3_ADV_X( bi,bj,k, dTtracerLev(k),
216         I           uTrans, uVel, maskW(1-Olx,1-Oly,k,bi,bj), localT,
217         O           af, myThid )
218             ELSE
219              CALL GAD_DST3FL_ADV_X( bi,bj,k, dTtracerLev(k),
220         I           uTrans, uVel, maskW(1-Olx,1-Oly,k,bi,bj), localT,
221         O           af, myThid )
222             ENDIF
223            ELSE
224             STOP 'GAD_CALC_RHS: Bad advectionScheme (X)'
225            ENDIF
226            DO j=1-Oly,sNy+Oly
227             DO i=1-Olx,sNx+Olx
228              fZon(i,j) = fZon(i,j) + af(i,j)
229             ENDDO
230            ENDDO
231    #ifdef ALLOW_DIAGNOSTICS
232            IF ( useDiagnostics ) THEN
233              diagName = 'ADVx'//diagSufx
234              CALL DIAGNOSTICS_FILL(af,diagName, k,1, 2,bi,bj, myThid)
235            ENDIF
236    #endif
237        ENDIF        ENDIF
238    
239  C-    Diffusive flux in X  C-    Diffusive flux in X
# Line 191  C-    Diffusive flux in X Line 247  C-    Diffusive flux in X
247         ENDDO         ENDDO
248        ENDIF        ENDIF
249    
250    C-    Add bi-harmonic diffusive flux in X
251          IF (diffK4 .NE. 0.) THEN
252           CALL GAD_BIHARM_X(bi,bj,k,xA,df4,diffK4,df,myThid)
253          ENDIF
254    
255  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
256  C-    GM/Redi flux in X  C-    GM/Redi flux in X
257        IF (useGMRedi) THEN        IF (useGMRedi) THEN
# Line 208  C *note* should update GMREDI_XTRANSPORT Line 269  C *note* should update GMREDI_XTRANSPORT
269         ENDDO         ENDDO
270        ENDDO        ENDDO
271    
272  C-    Bi-harmonic duffusive flux in X  #ifdef ALLOW_DIAGNOSTICS
273        IF (diffK4 .NE. 0.) THEN  C-    Diagnostics of Tracer flux in X dir (mainly Diffusive term),
274         CALL GAD_BIHARM_X(bi,bj,k,xA,df4,diffK4,df,myThid)  C       excluding advective terms:
275         DO j=1-Oly,sNy+Oly        IF ( useDiagnostics .AND.
276          DO i=1-Olx,sNx+Olx       &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. useGMRedi) ) THEN
277           fZon(i,j) = fZon(i,j) + df(i,j)            diagName = 'DIFx'//diagSufx
278          ENDDO            CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
        ENDDO  
279        ENDIF        ENDIF
280    #endif
281    
282  C--   Initialize net flux in Y direction  C--   Initialize net flux in Y direction
283        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
# Line 227  C--   Initialize net flux in Y direction Line 288  C--   Initialize net flux in Y direction
288    
289  C-    Advective flux in Y  C-    Advective flux in Y
290        IF (calcAdvection) THEN        IF (calcAdvection) THEN
291        IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN          IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
292         CALL GAD_C2_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)            CALL GAD_C2_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)
293        ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST
294         CALL GAD_FLUXLIMIT_ADV_Y(       &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
295       &       bi,bj,k,deltaTtracer,vTrans,vVel,localT,af,myThid)            CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme,
296        ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN       I            dTtracerLev(k), vTrans, vVel, localT,
297         CALL GAD_U3_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)       O            af, myThid )
298        ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
299         CALL GAD_C4_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)            CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, dTtracerLev(k),
300        ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN       I            vTrans, vVel, maskS(1-Olx,1-Oly,k,bi,bj), localT,
301         CALL GAD_DST3_ADV_Y(       O            af, myThid )
302       &       bi,bj,k,deltaTtracer,vTrans,vVel,localT,af,myThid)          ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
303        ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN            CALL GAD_U3_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)
304         CALL GAD_DST3FL_ADV_Y(          ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
305       &       bi,bj,k,deltaTtracer,vTrans,vVel,localT,af,myThid)            CALL GAD_C4_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)
306        ELSE          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
307         STOP 'GAD_CALC_RHS: Bad advectionScheme (Y)'            CALL GAD_DST3_ADV_Y( bi,bj,k, dTtracerLev(k),
308        ENDIF       I            vTrans, vVel, maskS(1-Olx,1-Oly,k,bi,bj), localT,
309        DO j=1-Oly,sNy+Oly       O            af, myThid )
310         DO i=1-Olx,sNx+Olx          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
311          fMer(i,j) = fMer(i,j) + af(i,j)           IF ( inAdMode ) THEN
312         ENDDO  cph This block is to trick the adjoint:
313        ENDDO  cph IF inAdExact=.FALSE., we want to use DST3
314    cph with limiters in forward, but without limiters in reverse.
315              CALL GAD_DST3_ADV_Y( bi,bj,k, dTtracerLev(k),
316         I           vTrans, vVel, maskS(1-Olx,1-Oly,k,bi,bj), localT,
317         O           af, myThid )
318             ELSE
319              CALL GAD_DST3FL_ADV_Y( bi,bj,k, dTtracerLev(k),
320         I           vTrans, vVel, maskS(1-Olx,1-Oly,k,bi,bj), localT,
321         O           af, myThid )
322             ENDIF
323            ELSE
324              STOP 'GAD_CALC_RHS: Bad advectionScheme (Y)'
325            ENDIF
326            DO j=1-Oly,sNy+Oly
327             DO i=1-Olx,sNx+Olx
328              fMer(i,j) = fMer(i,j) + af(i,j)
329             ENDDO
330            ENDDO
331    #ifdef ALLOW_DIAGNOSTICS
332            IF ( useDiagnostics ) THEN
333              diagName = 'ADVy'//diagSufx
334              CALL DIAGNOSTICS_FILL(af,diagName, k,1, 2,bi,bj, myThid)
335            ENDIF
336    #endif
337        ENDIF        ENDIF
338    
339  C-    Diffusive flux in Y  C-    Diffusive flux in Y
# Line 263  C-    Diffusive flux in Y Line 347  C-    Diffusive flux in Y
347         ENDDO         ENDDO
348        ENDIF        ENDIF
349    
350    C-    Add bi-harmonic flux in Y
351          IF (diffK4 .NE. 0.) THEN
352           CALL GAD_BIHARM_Y(bi,bj,k,yA,df4,diffK4,df,myThid)
353          ENDIF
354    
355  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
356  C-    GM/Redi flux in Y  C-    GM/Redi flux in Y
357        IF (useGMRedi) THEN        IF (useGMRedi) THEN
# Line 280  C *note* should update GMREDI_YTRANSPORT Line 369  C *note* should update GMREDI_YTRANSPORT
369         ENDDO         ENDDO
370        ENDDO        ENDDO
371    
372  C-    Bi-harmonic flux in Y  #ifdef ALLOW_DIAGNOSTICS
373        IF (diffK4 .NE. 0.) THEN  C-    Diagnostics of Tracer flux in Y dir (mainly Diffusive terms),
374         CALL GAD_BIHARM_Y(bi,bj,k,yA,df4,diffK4,df,myThid)  C       excluding advective terms:
375         DO j=1-Oly,sNy+Oly        IF ( useDiagnostics .AND.
376          DO i=1-Olx,sNx+Olx       &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. useGMRedi) ) THEN
377           fMer(i,j) = fMer(i,j) + df(i,j)            diagName = 'DIFy'//diagSufx
378          ENDDO            CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
        ENDDO  
       ENDIF  
   
 #ifdef NONLIN_FRSURF  
 C--   Compute vertical flux fVerT(kDown) at interface k+1 (between k & k+1):  
       IF ( calcAdvection .AND. K.EQ.Nr .AND.  
      &     useRealFreshWaterFlux .AND.  
      &     buoyancyRelation .EQ. 'OCEANICP' ) THEN    
        DO j=1-Oly,sNy+Oly  
         DO i=1-Olx,sNx+Olx  
          fVerT(i,j,kDown) = convertEmP2rUnit*PmEpR(i,j,bi,bj)  
      &     *rA(i,j,bi,bj)*maskC(i,j,k,bi,bj)*Tracer(i,j,k,bi,bj)  
         ENDDO  
        ENDDO  
379        ENDIF        ENDIF
380  #endif /* NONLIN_FRSURF */  #endif
381    
382  C--   Compute vertical flux fVerT(kUp) at interface k (between k-1 & k):  C--   Compute vertical flux fVerT(kUp) at interface k (between k-1 & k):
383  C-    Advective flux in R  C-    Advective flux in R
384        IF (calcAdvection) THEN  #ifdef ALLOW_AIM
385  C     Note: wVel needs to be masked  C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr
386        IF (K.GE.2) THEN        IF (calcAdvection .AND. .NOT.implicitAdvection .AND. K.GE.2 .AND.
387         &     (.NOT.useAIM .OR.tracerIdentity.NE.GAD_SALINITY .OR.K.LT.Nr)
388         &   ) THEN
389    #else
390          IF (calcAdvection .AND. .NOT.implicitAdvection .AND. K.GE.2) THEN
391    #endif
392  C-    Compute vertical advective flux in the interior:  C-    Compute vertical advective flux in the interior:
393         IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN          IF (vertAdvecScheme.EQ.ENUM_CENTERED_2ND) THEN
394          CALL GAD_C2_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)            CALL GAD_C2_ADV_R(bi,bj,k,rTrans,Tracer,af,myThid)
395         ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF ( vertAdvecScheme.EQ.ENUM_UPWIND_1RST
396          CALL GAD_FLUXLIMIT_ADV_R(       &          .OR. vertAdvecScheme.EQ.ENUM_DST2 ) THEN
397       &       bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)            CALL GAD_DST2U1_ADV_R( bi,bj,k, vertAdvecScheme,
398         ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN       I         dTtracerLev(k),rTrans,wVel,Tracer(1-Olx,1-Oly,1,bi,bj),
399          CALL GAD_U3_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)       O         af, myThid )
400         ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN          ELSEIF (vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN
401          CALL GAD_C4_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)            CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k,
402         ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN       I         dTtracerLev(k),rTrans,wVel,Tracer(1-Olx,1-Oly,1,bi,bj),
403          CALL GAD_DST3_ADV_R(       O         af, myThid )
404       &       bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)          ELSEIF (vertAdvecScheme.EQ.ENUM_UPWIND_3RD ) THEN
405         ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN            CALL GAD_U3_ADV_R(bi,bj,k,rTrans,Tracer,af,myThid)
406          CALL GAD_DST3FL_ADV_R(          ELSEIF (vertAdvecScheme.EQ.ENUM_CENTERED_4TH) THEN
407       &       bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)            CALL GAD_C4_ADV_R(bi,bj,k,rTrans,Tracer,af,myThid)
408         ELSE          ELSEIF (vertAdvecScheme.EQ.ENUM_DST3 ) THEN
409          STOP 'GAD_CALC_RHS: Bad advectionScheme (R)'            CALL GAD_DST3_ADV_R( bi,bj,k,
410         ENDIF       I         dTtracerLev(k),rTrans,wVel,Tracer(1-Olx,1-Oly,1,bi,bj),
411  C-    Surface "correction" term at k>1 :       O         af, myThid )
412         DO j=1-Oly,sNy+Oly          ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
413          DO i=1-Olx,sNx+Olx  cph This block is to trick the adjoint:
414           af(i,j) = af(i,j)  cph IF inAdExact=.FALSE., we want to use DST3
415       &           + (maskC(i,j,k,bi,bj)-maskC(i,j,k-1,bi,bj))*  cph with limiters in forward, but without limiters in reverse.
416       &             rTrans(i,j)*Tracer(i,j,k,bi,bj)            IF ( inAdMode ) THEN
417          ENDDO             CALL GAD_DST3_ADV_R( bi,bj,k,
418         ENDDO       I         dTtracerLev(k),rTrans,wVel,Tracer(1-Olx,1-Oly,1,bi,bj),
419        ELSE       O         af, myThid )
420  C-    Surface "correction" term at k=1 :            ELSE
421         DO j=1-Oly,sNy+Oly             CALL GAD_DST3FL_ADV_R( bi,bj,k,
422          DO i=1-Olx,sNx+Olx       I         dTtracerLev(k),rTrans,wVel,Tracer(1-Olx,1-Oly,1,bi,bj),
423           af(i,j) = rTrans(i,j)*Tracer(i,j,k,bi,bj)       O         af, myThid )
424          ENDDO            ENDIF
425         ENDDO          ELSE
426        ENDIF            STOP 'GAD_CALC_RHS: Bad vertAdvecScheme (R)'
427  C-    add the advective flux to fVerT          ENDIF
428        DO j=1-Oly,sNy+Oly  C-     add the advective flux to fVerT
429         DO i=1-Olx,sNx+Olx          DO j=1-Oly,sNy+Oly
430          fVerT(i,j,kUp) = fVerT(i,j,kUp) + af(i,j)           DO i=1-Olx,sNx+Olx
431         ENDDO            fVerT(i,j,kUp) = fVerT(i,j,kUp) + af(i,j)
432        ENDDO           ENDDO
433            ENDDO
434    #ifdef ALLOW_DIAGNOSTICS
435            IF ( useDiagnostics ) THEN
436              diagName = 'ADVr'//diagSufx
437              CALL DIAGNOSTICS_FILL(af,diagName, k,1, 2,bi,bj, myThid)
438    C- note: needs to explicitly increment the counter since DIAGNOSTICS_FILL
439    C        does it only if k=1 (never the case here)
440              IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT(diagName,bi,bj,myThid)
441            ENDIF
442    #endif
443        ENDIF        ENDIF
444    
445  C-    Diffusive flux in R  C-    Diffusive flux in R
# Line 362  C           boundary condition. Line 452  C           boundary condition.
452          ENDDO          ENDDO
453         ENDDO         ENDDO
454        ELSE        ELSE
455         CALL GAD_DIFF_R(bi,bj,k,KappaRT,tracer,df,myThid)         CALL GAD_DIFF_R(bi,bj,k,KappaR,Tracer,df,myThid)
456        ENDIF        ENDIF
457    
458  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
# Line 383  C *note* should update GMREDI_RTRANSPORT Line 473  C *note* should update GMREDI_RTRANSPORT
473         ENDDO         ENDDO
474        ENDDO        ENDDO
475    
476    #ifdef ALLOW_DIAGNOSTICS
477    C-    Diagnostics of Tracer flux in R dir (mainly Diffusive terms),
478    C       Explicit terms only & excluding advective terms:
479          IF ( useDiagnostics .AND.
480         &    (.NOT.implicitDiffusion .OR. useGMRedi) ) THEN
481              diagName = 'DFrE'//diagSufx
482              CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
483          ENDIF
484    #endif
485    
486  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
487  C-    Add non local KPP transport term (ghat) to diffusive T flux.  C-    Set non local KPP transport term (ghat):
488        IF (useKPP) THEN        IF ( useKPP .AND. k.GE.2 ) THEN
489         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
490          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
491           df(i,j) = 0. _d 0           df(i,j) = 0. _d 0
492          ENDDO          ENDDO
493         ENDDO         ENDDO
494         IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN         IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
 C *note* should update KPP_TRANSPORT_T to set df  *aja*  
495          CALL KPP_TRANSPORT_T(          CALL KPP_TRANSPORT_T(
496       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,
497       I     KappaRT,       O     df )
      U     df )  
498         ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN         ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
499          CALL KPP_TRANSPORT_S(          CALL KPP_TRANSPORT_S(
500       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,
501       I     KappaRT,       O     df )
502       U     df )  #ifdef ALLOW_PTRACERS
503           ELSEIF (tracerIdentity .GE. GAD_TR1) THEN
504            CALL KPP_TRANSPORT_PTR(
505         I     iMin,iMax,jMin,jMax,bi,bj,k,km1,
506         I     tracerIdentity-GAD_TR1+1,
507         O     df )
508    #endif
509         ELSE         ELSE
510            PRINT*,'invalid tracer indentity: ', tracerIdentity
511          STOP 'GAD_CALC_RHS: Ooops'          STOP 'GAD_CALC_RHS: Ooops'
512         ENDIF         ENDIF
513         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
# Line 417  C--   Divergence of fluxes Line 522  C--   Divergence of fluxes
522        DO j=1-Oly,sNy+Oly-1        DO j=1-Oly,sNy+Oly-1
523         DO i=1-Olx,sNx+Olx-1         DO i=1-Olx,sNx+Olx-1
524          gTracer(i,j,k,bi,bj)=gTracer(i,j,k,bi,bj)          gTracer(i,j,k,bi,bj)=gTracer(i,j,k,bi,bj)
525       &   -_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)
526       &    *recip_rA(i,j,bi,bj)       &   *( (fZon(i+1,j)-fZon(i,j))
527       &    *(       &     +(fMer(i,j+1)-fMer(i,j))
528       &    +( fZon(i+1,j)-fZon(i,j) )       &     +(fVerT(i,j,kDown)-fVerT(i,j,kUp))*rkSign
529       &    +( fMer(i,j+1)-fMer(i,j) )       &     -localT(i,j)*( (uTrans(i+1,j)-uTrans(i,j))
530       &    +( fVerT(i,j,kUp)-fVerT(i,j,kDown) )*rkFac       &                   +(vTrans(i,j+1)-vTrans(i,j))
531         &                   +(rTransKp1(i,j)-rTrans(i,j))*rAdvFac
532         &                  )*advFac
533       &    )       &    )
534         ENDDO         ENDDO
535        ENDDO        ENDDO
536    
537  #ifdef NONLIN_FRSURF  #ifdef ALLOW_DEBUG
538  C-- account for 3.D divergence of the flow in rStar coordinate:        IF ( debugLevel .GE. debLevB
539        IF (calcAdvection .AND. select_rStar.GT.0) THEN       &   .AND. tracerIdentity.EQ.GAD_TEMPERATURE
540         DO j=1-Oly,sNy+Oly-1       &   .AND. k.EQ.2 .AND. myIter.EQ.1+nIter0
541          DO i=1-Olx,sNx+Olx-1       &   .AND. nPx.EQ.1 .AND. nPy.EQ.1
542           gTracer(i,j,k,bi,bj) = gTracer(i,j,k,bi,bj)       &   .AND. useCubedSphereExchange ) THEN
543       &     - (rStarExpC(i,j,bi,bj) - 1. _d 0)/deltaTfreesurf          CALL DEBUG_CS_CORNER_UV( ' fZon,fMer from GAD_CALC_RHS',
544       &       *tracer(i,j,k,bi,bj)*maskC(i,j,k,bi,bj)       &             fZon,fMer, k, standardMessageUnit,bi,bj,myThid )
         ENDDO  
        ENDDO  
545        ENDIF        ENDIF
546        IF (calcAdvection .AND. select_rStar.LT.0) THEN  #endif /* ALLOW_DEBUG */
547         DO j=1-Oly,sNy+Oly-1  
         DO i=1-Olx,sNx+Olx-1  
          gTracer(i,j,k,bi,bj) = gTracer(i,j,k,bi,bj)  
      &     - rStarDhCDt(i,j,bi,bj)  
      &       *tracer(i,j,k,bi,bj)*maskC(i,j,k,bi,bj)  
         ENDDO  
        ENDDO  
       ENDIF  
 #endif /* NONLIN_FRSURF */  
         
   
548        RETURN        RETURN
549        END        END
Legend:
Removed from v.1.17  
changed lines
  Added in v.1.37
  ViewVC Help
Powered by ViewVC 1.1.22