/[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.19 by dimitri, Mon Aug 4 22:53:42 2003 UTC revision 1.47 by jmc, Thu May 3 21:34:39 2007 UTC
# Line 7  CBOP Line 7  CBOP
7  C !ROUTINE: GAD_CALC_RHS  C !ROUTINE: GAD_CALC_RHS
8    
9  C !INTERFACE: ==========================================================  C !INTERFACE: ==========================================================
10        SUBROUTINE GAD_CALC_RHS(        SUBROUTINE GAD_CALC_RHS(
11       I           bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,       I           bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,
12       I           xA,yA,uTrans,vTrans,rTrans,maskUp,       I           xA, yA, maskUp, uFld, vFld, wFld,
13       I           diffKh, diffK4, KappaRT, Tracer,       I           uTrans, vTrans, rTrans, rTransKp1,
14       I           tracerIdentity, advectionScheme, calcAdvection,       I           diffKh, diffK4, KappaR, TracerN, TracAB,
15         I           tracerIdentity, advectionScheme, vertAdvecScheme,
16         I           calcAdvection, implicitAdvection, applyAB_onTracer,
17       U           fVerT, gTracer,       U           fVerT, gTracer,
18       I           myThid )       I           myTime, myIter, myThid )
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"
 #include "DYNVARS.h"  
45  #include "SURFACE.h"  #include "SURFACE.h"
46  #include "GAD.h"  #include "GAD.h"
 #ifdef ALLOW_PTRACERS  
 #include "PTRACERS_OPTIONS.h"  
 #include "PTRACERS.h"  
 #endif  
47    
48  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
49  #include "tamc.h"  #include "tamc.h"
# Line 54  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 k                :: vertical index
58  C  kp1                  :: =k+1 for k<Nr, =Nr for k=Nr  C kM1              :: =k-1 for k>1, =1 for k=1
59  C  xA,yA                :: areas of X and Y face of tracer cells  C kUp              :: index into 2 1/2D array, toggles between 1|2
60  C  uTrans,vTrans,rTrans :: 2-D arrays of volume transports at U,V and W points  C kDown            :: index into 2 1/2D array, toggles between 2|1
61  C  maskUp               :: 2-D array for mask at W points  C xA,yA            :: areas of X and Y face of tracer cells
62  C  diffKh               :: horizontal diffusion coefficient  C maskUp           :: 2-D array for mask at W points
63  C  diffK4               :: bi-harmonic diffusion coefficient  C uFld,vFld,wFld   :: Local copy of velocity field (3 components)
64  C  KappaRT              :: 3-D array for vertical diffusion coefficient  C uTrans,vTrans    :: 2-D arrays of volume transports at U,V points
65  C  Tracer               :: tracer field  C rTrans           :: 2-D arrays of volume transports at W points
66  C  tracerIdentity       :: identifier for the tracer (required for KPP and GM)  C rTransKp1        :: 2-D array of volume trans at W pts, interf k+1
67  C  advectionScheme      :: advection scheme to use  C diffKh           :: horizontal diffusion coefficient
68  C  calcAdvection        :: =False if Advec terms computed with multiDim scheme  C diffK4           :: bi-harmonic diffusion coefficient
69  C  myThid               :: thread number  C KappaR           :: 2-D array for vertical diffusion coefficient, interf k
70    C TracerN          :: tracer field @ time-step n (Note: only used
71    C                     if applying AB on tracer field rather than on tendency gTr)
72    C TracAB           :: current tracer field (@ time-step n if applying AB on gTr
73    C                     or extrapolated fwd in time to n+1/2 if applying AB on Tr)
74    C tracerIdentity   :: tracer identifier (required for KPP,GM)
75    C advectionScheme  :: advection scheme to use (Horizontal plane)
76    C vertAdvecScheme  :: advection scheme to use (Vertical direction)
77    C calcAdvection    :: =False if Advec computed with multiDim scheme
78    C implicitAdvection:: =True if vertical Advec computed implicitly
79    C applyAB_onTracer :: apply Adams-Bashforth on Tracer (rather than on gTr)
80    C myTime           :: current time
81    C myIter           :: iteration number
82    C myThid           :: thread number
83        INTEGER bi,bj,iMin,iMax,jMin,jMax        INTEGER bi,bj,iMin,iMax,jMin,jMax
84        INTEGER k,kUp,kDown,kM1        INTEGER k,kUp,kDown,kM1
85        _RS xA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS xA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
86        _RS yA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS yA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
87          _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
88          _RL uFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
89          _RL vFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
90          _RL wFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
91        _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
92        _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
93        _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
94        _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
95        _RL diffKh, diffK4        _RL diffKh, diffK4
96        _RL KappaRT(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL KappaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
97        _RL Tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL TracerN(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
98          _RL TracAB (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
99        INTEGER tracerIdentity        INTEGER tracerIdentity
100        INTEGER advectionScheme        INTEGER advectionScheme, vertAdvecScheme
101        LOGICAL calcAdvection        LOGICAL calcAdvection
102        INTEGER myThid        LOGICAL implicitAdvection, applyAB_onTracer
103          _RL     myTime
104          INTEGER myIter, myThid
105    
106  C !OUTPUT PARAMETERS: ==================================================  C !OUTPUT PARAMETERS: ==================================================
107  C  gTracer              :: tendancy array  C gTracer          :: tendency array
108  C  fVerT                :: 2 1/2D arrays for vertical advective flux  C fVerT            :: 2 1/2D arrays for vertical advective flux
109        _RL gTracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL gTracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
110        _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
111    
112  C !LOCAL VARIABLES: ====================================================  C !LOCAL VARIABLES: ====================================================
113  C  i,j                  :: loop indices  C i,j              :: loop indices
114  C  df4                  :: used for storing del^2 T for bi-harmonic term  C df4              :: used for storing del^2 T for bi-harmonic term
115  C  fZon                 :: zonal flux  C fZon             :: zonal flux
116  C  fmer                 :: meridional flux  C fMer             :: meridional flux
117  C  af                   :: advective flux  C af               :: advective flux
118  C  df                   :: diffusive flux  C df               :: diffusive flux
119  C  localT               :: local copy of tracer field  C localT           :: local copy of tracer field
120    C locABT           :: local copy of (AB-extrapolated) tracer field
121    #ifdef ALLOW_DIAGNOSTICS
122          CHARACTER*8 diagName
123          CHARACTER*4 GAD_DIAG_SUFX, diagSufx
124          EXTERNAL    GAD_DIAG_SUFX
125    #endif
126        INTEGER i,j        INTEGER i,j
127        _RL df4   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL df4   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
128        _RL fZon  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fZon  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 107  C  localT               :: local copy of Line 130  C  localT               :: local copy of
130        _RL af    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL af    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
131        _RL df    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL df    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
132        _RL localT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL localT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
133          _RL locABT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
134          _RL advFac, rAdvFac
135  CEOP  CEOP
136    
137  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
# Line 115  C--   the kDown is still required Line 140  C--   the kDown is still required
140        fVerT(1,1,kDown) = fVerT(1,1,kDown)        fVerT(1,1,kDown) = fVerT(1,1,kDown)
141  #endif  #endif
142    
143    #ifdef ALLOW_DIAGNOSTICS
144    C--   Set diagnostic suffix for the current tracer
145          IF ( useDiagnostics ) THEN
146            diagSufx = GAD_DIAG_SUFX( tracerIdentity, myThid )
147          ENDIF
148    #endif
149    
150          advFac  = 0. _d 0
151          IF (calcAdvection) advFac = 1. _d 0
152          rAdvFac = rkSign*advFac
153          IF (implicitAdvection) rAdvFac = 0. _d 0
154    
155        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
156         DO i=1-OLx,sNx+OLx         DO i=1-OLx,sNx+OLx
157          fZon(i,j)      = 0. _d 0          fZon(i,j)      = 0. _d 0
# Line 122  C--   the kDown is still required Line 159  C--   the kDown is still required
159          fVerT(i,j,kUp) = 0. _d 0          fVerT(i,j,kUp) = 0. _d 0
160          df(i,j)        = 0. _d 0          df(i,j)        = 0. _d 0
161          df4(i,j)       = 0. _d 0          df4(i,j)       = 0. _d 0
         localT(i,j)    = 0. _d 0  
162         ENDDO         ENDDO
163        ENDDO        ENDDO
164    
165  C--   Make local copy of tracer array  C--   Make local copy of tracer array
166        DO j=1-OLy,sNy+OLy        IF ( applyAB_onTracer ) THEN
167         DO i=1-OLx,sNx+OLx          DO j=1-OLy,sNy+OLy
168          localT(i,j)=tracer(i,j,k,bi,bj)           DO i=1-OLx,sNx+OLx
169         ENDDO            localT(i,j)=TracerN(i,j,k,bi,bj)
170        ENDDO            locABT(i,j)= TracAB(i,j,k,bi,bj)
171             ENDDO
172            ENDDO
173          ELSE
174            DO j=1-OLy,sNy+OLy
175             DO i=1-OLx,sNx+OLx
176              localT(i,j)= TracAB(i,j,k,bi,bj)
177              locABT(i,j)= TracAB(i,j,k,bi,bj)
178             ENDDO
179            ENDDO
180          ENDIF
181    
182  C--   Unless we have already calculated the advection terms we initialize  C--   Unless we have already calculated the advection terms we initialize
183  C     the tendency to zero.  C     the tendency to zero.
184        IF (calcAdvection) THEN  C     <== now done earlier at the beginning of thermodynamics.
185         DO j=1-Oly,sNy+Oly  c     IF (calcAdvection) THEN
186          DO i=1-Olx,sNx+Olx  c      DO j=1-Oly,sNy+Oly
187           gTracer(i,j,k,bi,bj)=0. _d 0  c       DO i=1-Olx,sNx+Olx
188          ENDDO  c        gTracer(i,j,k,bi,bj)=0. _d 0
189         ENDDO  c       ENDDO
190        ENDIF  c      ENDDO
191    c     ENDIF
192    
193  C--   Pre-calculate del^2 T if bi-harmonic coefficient is non-zero  C--   Pre-calculate del^2 T if bi-harmonic coefficient is non-zero
194        IF (diffK4 .NE. 0.) THEN        IF (diffK4 .NE. 0.) THEN
# Line 159  C--   Initialize net flux in X direction Line 206  C--   Initialize net flux in X direction
206    
207  C-    Advective flux in X  C-    Advective flux in X
208        IF (calcAdvection) THEN        IF (calcAdvection) THEN
209        IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN          IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
210         CALL GAD_C2_ADV_X(bi,bj,k,uTrans,localT,af,myThid)            CALL GAD_C2_ADV_X(bi,bj,k,uTrans,locABT,af,myThid)
211        ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST
212         CALL GAD_FLUXLIMIT_ADV_X(       &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
213       &      bi,bj,k,deltaTtracer,uTrans,uVel,localT,af,myThid)            CALL GAD_DST2U1_ADV_X( bi,bj,k, advectionScheme, .TRUE.,
214        ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN       I            dTtracerLev(k), uTrans, uFld, locABT,
215         CALL GAD_U3_ADV_X(bi,bj,k,uTrans,localT,af,myThid)       O            af, myThid )
216        ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
217         CALL GAD_C4_ADV_X(bi,bj,k,uTrans,localT,af,myThid)            CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, .TRUE., dTtracerLev(k),
218        ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN       I            uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,
219         CALL GAD_DST3_ADV_X(       O            af, myThid )
220       &       bi,bj,k,deltaTtracer,uTrans,uVel,localT,af,myThid)          ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
221        ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN            CALL GAD_U3_ADV_X(bi,bj,k,uTrans,locABT,af,myThid)
222         CALL GAD_DST3FL_ADV_X(          ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
223       &       bi,bj,k,deltaTtracer,uTrans,uVel,localT,af,myThid)            CALL GAD_C4_ADV_X(bi,bj,k,uTrans,locABT,af,myThid)
224        ELSE          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
225         STOP 'GAD_CALC_RHS: Bad advectionScheme (X)'            CALL GAD_DST3_ADV_X( bi,bj,k, .TRUE., dTtracerLev(k),
226        ENDIF       I            uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,
227        DO j=1-Oly,sNy+Oly       O            af, myThid )
228         DO i=1-Olx,sNx+Olx          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
229          fZon(i,j) = fZon(i,j) + af(i,j)           IF ( inAdMode ) THEN
230         ENDDO  cph This block is to trick the adjoint:
231        ENDDO  cph IF inAdExact=.FALSE., we want to use DST3
232    cph with limiters in forward, but without limiters in reverse.
233              CALL GAD_DST3_ADV_X( bi,bj,k, .TRUE., dTtracerLev(k),
234         I           uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,
235         O           af, myThid )
236             ELSE
237              CALL GAD_DST3FL_ADV_X( bi,bj,k, .TRUE., dTtracerLev(k),
238         I           uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,
239         O           af, myThid )
240             ENDIF
241    #ifndef ALLOW_AUTODIFF_TAMC
242            ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN
243              CALL GAD_OS7MP_ADV_X( bi,bj,k, .TRUE., dTtracerLev(k),
244         I            uTrans, uFld, maskW(1-Olx,1-Oly,k,bi,bj), locABT,
245         O            af, myThid )
246    #endif
247            ELSE
248             STOP 'GAD_CALC_RHS: Bad advectionScheme (X)'
249            ENDIF
250            DO j=1-Oly,sNy+Oly
251             DO i=1-Olx,sNx+Olx
252              fZon(i,j) = fZon(i,j) + af(i,j)
253             ENDDO
254            ENDDO
255    #ifdef ALLOW_DIAGNOSTICS
256            IF ( useDiagnostics ) THEN
257              diagName = 'ADVx'//diagSufx
258              CALL DIAGNOSTICS_FILL(af,diagName, k,1, 2,bi,bj, myThid)
259            ENDIF
260    #endif
261        ENDIF        ENDIF
262    
263  C-    Diffusive flux in X  C-    Diffusive flux in X
# Line 195  C-    Diffusive flux in X Line 271  C-    Diffusive flux in X
271         ENDDO         ENDDO
272        ENDIF        ENDIF
273    
274    C-    Add bi-harmonic diffusive flux in X
275          IF (diffK4 .NE. 0.) THEN
276           CALL GAD_BIHARM_X(bi,bj,k,xA,df4,diffK4,df,myThid)
277          ENDIF
278    
279  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
280  C-    GM/Redi flux in X  C-    GM/Redi flux in X
281        IF (useGMRedi) THEN        IF (useGMRedi) THEN
282  C *note* should update GMREDI_XTRANSPORT to use localT and set df  *aja*  C *note* should update GMREDI_XTRANSPORT to set df  *aja*
283          CALL GMREDI_XTRANSPORT(          IF ( applyAB_onTracer ) THEN
284       I     iMin,iMax,jMin,jMax,bi,bj,K,            CALL GMREDI_XTRANSPORT(
285       I     xA,Tracer,tracerIdentity,       I         iMin,iMax,jMin,jMax,bi,bj,k,
286       U     df,       I         xA,TracerN,tracerIdentity,
287       I     myThid)       U         df,
288         I         myThid)
289            ELSE
290              CALL GMREDI_XTRANSPORT(
291         I         iMin,iMax,jMin,jMax,bi,bj,k,
292         I         xA,TracAB, tracerIdentity,
293         U         df,
294         I         myThid)
295            ENDIF
296        ENDIF        ENDIF
297  #endif  #endif
298    C     anelastic: advect.fluxes are scaled by rhoFac but hor.diff. flx are not
299        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
300         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
301          fZon(i,j) = fZon(i,j) + df(i,j)          fZon(i,j) = fZon(i,j) + df(i,j)*rhoFacC(k)
302         ENDDO         ENDDO
303        ENDDO        ENDDO
304    
305  C-    Bi-harmonic duffusive flux in X  #ifdef ALLOW_DIAGNOSTICS
306        IF (diffK4 .NE. 0.) THEN  C-    Diagnostics of Tracer flux in X dir (mainly Diffusive term),
307         CALL GAD_BIHARM_X(bi,bj,k,xA,df4,diffK4,df,myThid)  C       excluding advective terms:
308         DO j=1-Oly,sNy+Oly        IF ( useDiagnostics .AND.
309          DO i=1-Olx,sNx+Olx       &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. useGMRedi) ) THEN
310           fZon(i,j) = fZon(i,j) + df(i,j)            diagName = 'DFxE'//diagSufx
311          ENDDO            CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
        ENDDO  
312        ENDIF        ENDIF
313    #endif
314    
315  C--   Initialize net flux in Y direction  C--   Initialize net flux in Y direction
316        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
# Line 231  C--   Initialize net flux in Y direction Line 321  C--   Initialize net flux in Y direction
321    
322  C-    Advective flux in Y  C-    Advective flux in Y
323        IF (calcAdvection) THEN        IF (calcAdvection) THEN
324        IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN          IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
325         CALL GAD_C2_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)            CALL GAD_C2_ADV_Y(bi,bj,k,vTrans,locABT,af,myThid)
326        ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST
327         CALL GAD_FLUXLIMIT_ADV_Y(       &          .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
328       &       bi,bj,k,deltaTtracer,vTrans,vVel,localT,af,myThid)            CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme, .TRUE.,
329        ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN       I            dTtracerLev(k), vTrans, vFld, locABT,
330         CALL GAD_U3_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)       O            af, myThid )
331        ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN          ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
332         CALL GAD_C4_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)            CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, .TRUE., dTtracerLev(k),
333        ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN       I            vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,
334         CALL GAD_DST3_ADV_Y(       O            af, myThid )
335       &       bi,bj,k,deltaTtracer,vTrans,vVel,localT,af,myThid)          ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN
336        ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN            CALL GAD_U3_ADV_Y(bi,bj,k,vTrans,locABT,af,myThid)
337         CALL GAD_DST3FL_ADV_Y(          ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
338       &       bi,bj,k,deltaTtracer,vTrans,vVel,localT,af,myThid)            CALL GAD_C4_ADV_Y(bi,bj,k,vTrans,locABT,af,myThid)
339        ELSE          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
340         STOP 'GAD_CALC_RHS: Bad advectionScheme (Y)'            CALL GAD_DST3_ADV_Y( bi,bj,k, .TRUE., dTtracerLev(k),
341        ENDIF       I            vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,
342        DO j=1-Oly,sNy+Oly       O            af, myThid )
343         DO i=1-Olx,sNx+Olx          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
344          fMer(i,j) = fMer(i,j) + af(i,j)           IF ( inAdMode ) THEN
345         ENDDO  cph This block is to trick the adjoint:
346        ENDDO  cph IF inAdExact=.FALSE., we want to use DST3
347    cph with limiters in forward, but without limiters in reverse.
348              CALL GAD_DST3_ADV_Y( bi,bj,k, .TRUE., dTtracerLev(k),
349         I           vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,
350         O           af, myThid )
351             ELSE
352              CALL GAD_DST3FL_ADV_Y( bi,bj,k, .TRUE., dTtracerLev(k),
353         I           vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,
354         O           af, myThid )
355             ENDIF
356    #ifndef ALLOW_AUTODIFF_TAMC
357            ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN
358              CALL GAD_OS7MP_ADV_Y( bi,bj,k, .TRUE., dTtracerLev(k),
359         I            vTrans, vFld, maskS(1-Olx,1-Oly,k,bi,bj), locABT,
360         O            af, myThid )
361    #endif
362            ELSE
363              STOP 'GAD_CALC_RHS: Bad advectionScheme (Y)'
364            ENDIF
365            DO j=1-Oly,sNy+Oly
366             DO i=1-Olx,sNx+Olx
367              fMer(i,j) = fMer(i,j) + af(i,j)
368             ENDDO
369            ENDDO
370    #ifdef ALLOW_DIAGNOSTICS
371            IF ( useDiagnostics ) THEN
372              diagName = 'ADVy'//diagSufx
373              CALL DIAGNOSTICS_FILL(af,diagName, k,1, 2,bi,bj, myThid)
374            ENDIF
375    #endif
376        ENDIF        ENDIF
377    
378  C-    Diffusive flux in Y  C-    Diffusive flux in Y
# Line 267  C-    Diffusive flux in Y Line 386  C-    Diffusive flux in Y
386         ENDDO         ENDDO
387        ENDIF        ENDIF
388    
389    C-    Add bi-harmonic flux in Y
390          IF (diffK4 .NE. 0.) THEN
391           CALL GAD_BIHARM_Y(bi,bj,k,yA,df4,diffK4,df,myThid)
392          ENDIF
393    
394  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
395  C-    GM/Redi flux in Y  C-    GM/Redi flux in Y
396        IF (useGMRedi) THEN        IF (useGMRedi) THEN
397  C *note* should update GMREDI_YTRANSPORT to use localT and set df  *aja*  C *note* should update GMREDI_YTRANSPORT to set df  *aja*
398         CALL GMREDI_YTRANSPORT(          IF ( applyAB_onTracer ) THEN
399       I     iMin,iMax,jMin,jMax,bi,bj,K,            CALL GMREDI_YTRANSPORT(
400       I     yA,Tracer,tracerIdentity,       I         iMin,iMax,jMin,jMax,bi,bj,k,
401       U     df,       I         yA,TracerN,tracerIdentity,
402       I     myThid)       U         df,
403         I         myThid)
404            ELSE
405              CALL GMREDI_YTRANSPORT(
406         I         iMin,iMax,jMin,jMax,bi,bj,k,
407         I         yA,TracAB, tracerIdentity,
408         U         df,
409         I         myThid)
410            ENDIF
411        ENDIF        ENDIF
412  #endif  #endif
413    C     anelastic: advect.fluxes are scaled by rhoFac but hor.diff. flx are not
414        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
415         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
416          fMer(i,j) = fMer(i,j) + df(i,j)          fMer(i,j) = fMer(i,j) + df(i,j)*rhoFacC(k)
417         ENDDO         ENDDO
418        ENDDO        ENDDO
419    
420  C-    Bi-harmonic flux in Y  #ifdef ALLOW_DIAGNOSTICS
421        IF (diffK4 .NE. 0.) THEN  C-    Diagnostics of Tracer flux in Y dir (mainly Diffusive terms),
422         CALL GAD_BIHARM_Y(bi,bj,k,yA,df4,diffK4,df,myThid)  C       excluding advective terms:
423         DO j=1-Oly,sNy+Oly        IF ( useDiagnostics .AND.
424          DO i=1-Olx,sNx+Olx       &    (diffKh.NE.0. .OR. diffK4 .NE.0. .OR. useGMRedi) ) THEN
425           fMer(i,j) = fMer(i,j) + df(i,j)            diagName = 'DFyE'//diagSufx
426          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  
427        ENDIF        ENDIF
428  #endif /* NONLIN_FRSURF */  #endif
429    
430  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):
431  C-    Advective flux in R  C-    Advective flux in R
432        IF (calcAdvection) THEN  #ifdef ALLOW_AIM
433  C     Note: wVel needs to be masked  C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr
434        IF (K.GE.2) THEN        IF (calcAdvection .AND. .NOT.implicitAdvection .AND. k.GE.2 .AND.
435         &     (.NOT.useAIM .OR.tracerIdentity.NE.GAD_SALINITY .OR.k.LT.Nr)
436         &   ) THEN
437    #else
438          IF (calcAdvection .AND. .NOT.implicitAdvection .AND. k.GE.2) THEN
439    #endif
440  C-    Compute vertical advective flux in the interior:  C-    Compute vertical advective flux in the interior:
441         IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN          IF (vertAdvecScheme.EQ.ENUM_CENTERED_2ND) THEN
442          CALL GAD_C2_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)            CALL GAD_C2_ADV_R(bi,bj,k,rTrans,TracAB,af,myThid)
443         ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN          ELSEIF ( vertAdvecScheme.EQ.ENUM_UPWIND_1RST
444          CALL GAD_FLUXLIMIT_ADV_R(       &          .OR. vertAdvecScheme.EQ.ENUM_DST2 ) THEN
445       &       bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)            CALL GAD_DST2U1_ADV_R( bi,bj,k, vertAdvecScheme,
446         ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD ) THEN       I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),
447          CALL GAD_U3_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)       O         af, myThid )
448         ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN          ELSEIF (vertAdvecScheme.EQ.ENUM_FLUX_LIMIT) THEN
449          CALL GAD_C4_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)            CALL GAD_FLUXLIMIT_ADV_R( bi,bj,k,
450         ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN       I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),
451          CALL GAD_DST3_ADV_R(       O         af, myThid )
452       &       bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)          ELSEIF (vertAdvecScheme.EQ.ENUM_UPWIND_3RD ) THEN
453         ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN            CALL GAD_U3_ADV_R(bi,bj,k,rTrans,TracAB,af,myThid)
454          CALL GAD_DST3FL_ADV_R(          ELSEIF (vertAdvecScheme.EQ.ENUM_CENTERED_4TH) THEN
455       &       bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)            CALL GAD_C4_ADV_R(bi,bj,k,rTrans,TracAB,af,myThid)
456         ELSE          ELSEIF (vertAdvecScheme.EQ.ENUM_DST3 ) THEN
457          STOP 'GAD_CALC_RHS: Bad advectionScheme (R)'            CALL GAD_DST3_ADV_R( bi,bj,k,
458         ENDIF       I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),
459  C-    Surface "correction" term at k>1 :       O         af, myThid )
460         DO j=1-Oly,sNy+Oly          ELSEIF (vertAdvecScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
461          DO i=1-Olx,sNx+Olx  cph This block is to trick the adjoint:
462           af(i,j) = af(i,j)  cph IF inAdExact=.FALSE., we want to use DST3
463       &           + (maskC(i,j,k,bi,bj)-maskC(i,j,k-1,bi,bj))*  cph with limiters in forward, but without limiters in reverse.
464       &             rTrans(i,j)*Tracer(i,j,k,bi,bj)            IF ( inAdMode ) THEN
465          ENDDO             CALL GAD_DST3_ADV_R( bi,bj,k,
466         ENDDO       I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),
467        ELSE       O         af, myThid )
468  C-    Surface "correction" term at k=1 :            ELSE
469         DO j=1-Oly,sNy+Oly             CALL GAD_DST3FL_ADV_R( bi,bj,k,
470          DO i=1-Olx,sNx+Olx       I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),
471           af(i,j) = rTrans(i,j)*Tracer(i,j,k,bi,bj)       O         af, myThid )
472          ENDDO            ENDIF
473         ENDDO  #ifndef ALLOW_AUTODIFF_TAMC
474        ENDIF          ELSEIF (vertAdvecScheme.EQ.ENUM_OS7MP ) THEN
475  C-    add the advective flux to fVerT             CALL GAD_OS7MP_ADV_R( bi,bj,k,
476        DO j=1-Oly,sNy+Oly       I         dTtracerLev(k),rTrans,wFld,TracAB(1-Olx,1-Oly,1,bi,bj),
477         DO i=1-Olx,sNx+Olx       O         af, myThid )
478          fVerT(i,j,kUp) = fVerT(i,j,kUp) + af(i,j)  #endif
479         ENDDO          ELSE
480        ENDDO            STOP 'GAD_CALC_RHS: Bad vertAdvecScheme (R)'
481            ENDIF
482    C-     add the advective flux to fVerT
483            DO j=1-Oly,sNy+Oly
484             DO i=1-Olx,sNx+Olx
485              fVerT(i,j,kUp) = fVerT(i,j,kUp) + af(i,j)
486             ENDDO
487            ENDDO
488    #ifdef ALLOW_DIAGNOSTICS
489            IF ( useDiagnostics ) THEN
490              diagName = 'ADVr'//diagSufx
491              CALL DIAGNOSTICS_FILL(af,diagName, k,1, 2,bi,bj, myThid)
492    C- note: needs to explicitly increment the counter since DIAGNOSTICS_FILL
493    C        does it only if k=1 (never the case here)
494              IF ( k.EQ.2 ) CALL DIAGNOSTICS_COUNT(diagName,bi,bj,myThid)
495            ENDIF
496    #endif
497        ENDIF        ENDIF
498    
499  C-    Diffusive flux in R  C-    Diffusive flux in R
# Line 366  C           boundary condition. Line 506  C           boundary condition.
506          ENDDO          ENDDO
507         ENDDO         ENDDO
508        ELSE        ELSE
509         CALL GAD_DIFF_R(bi,bj,k,KappaRT,tracer,df,myThid)         IF ( applyAB_onTracer ) THEN
510             CALL GAD_DIFF_R(bi,bj,k,KappaR,TracerN,df,myThid)
511           ELSE
512             CALL GAD_DIFF_R(bi,bj,k,KappaR,TracAB, df,myThid)
513           ENDIF
514        ENDIF        ENDIF
515    
516  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
517  C-    GM/Redi flux in R  C-    GM/Redi flux in R
518        IF (useGMRedi) THEN        IF (useGMRedi) THEN
519  C *note* should update GMREDI_RTRANSPORT to set df  *aja*  C *note* should update GMREDI_RTRANSPORT to set df  *aja*
520         CALL GMREDI_RTRANSPORT(          IF ( applyAB_onTracer ) THEN
521       I     iMin,iMax,jMin,jMax,bi,bj,K,            CALL GMREDI_RTRANSPORT(
522       I     Tracer,tracerIdentity,       I         iMin,iMax,jMin,jMax,bi,bj,k,
523       U     df,       I         TracerN,tracerIdentity,
524       I     myThid)       U         df,
525         I         myThid)
526            ELSE
527              CALL GMREDI_RTRANSPORT(
528         I         iMin,iMax,jMin,jMax,bi,bj,k,
529         I         TracAB, tracerIdentity,
530         U         df,
531         I         myThid)
532            ENDIF
533        ENDIF        ENDIF
534  #endif  #endif
535    
# Line 387  C *note* should update GMREDI_RTRANSPORT Line 539  C *note* should update GMREDI_RTRANSPORT
539         ENDDO         ENDDO
540        ENDDO        ENDDO
541    
542    #ifdef ALLOW_DIAGNOSTICS
543    C-    Diagnostics of Tracer flux in R dir (mainly Diffusive terms),
544    C       Explicit terms only & excluding advective terms:
545          IF ( useDiagnostics .AND.
546         &    (.NOT.implicitDiffusion .OR. useGMRedi) ) THEN
547              diagName = 'DFrE'//diagSufx
548              CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
549          ENDIF
550    #endif
551    
552  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
553  C-    Add non local KPP transport term (ghat) to diffusive T flux.  C-    Set non local KPP transport term (ghat):
554        IF (useKPP) THEN        IF ( useKPP .AND. k.GE.2 ) THEN
555         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
556          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
557           df(i,j) = 0. _d 0           df(i,j) = 0. _d 0
558          ENDDO          ENDDO
559         ENDDO         ENDDO
560         IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN         IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
 C *note* should update KPP_TRANSPORT_T to set df  *aja*  
561          CALL KPP_TRANSPORT_T(          CALL KPP_TRANSPORT_T(
562       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,       I           iMin,iMax,jMin,jMax,bi,bj,k,km1,
563       I     KappaRT,       O           df,
564       U     df )       I           myTime, myIter, myThid )
565         ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN         ELSEIF (tracerIdentity.EQ.GAD_SALINITY) THEN
566          CALL KPP_TRANSPORT_S(          CALL KPP_TRANSPORT_S(
567       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,       I           iMin,iMax,jMin,jMax,bi,bj,k,km1,
568       I     KappaRT,       O           df,
569       U     df )       I           myTime, myIter, myThid )
570  #ifdef ALLOW_PTRACERS  #ifdef ALLOW_PTRACERS
571         ELSEIF (tracerIdentity .GE. GAD_TR1 .AND.         ELSEIF (tracerIdentity .GE. GAD_TR1) THEN
      &         tracerIdentity .LE. (GAD_TR1+PTRACERS_numInUse-1)) THEN  
572          CALL KPP_TRANSPORT_PTR(          CALL KPP_TRANSPORT_PTR(
573       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,       I           iMin,iMax,jMin,jMax,bi,bj,k,km1,
574       I     tracerIdentity,KappaRT,       I           tracerIdentity-GAD_TR1+1,
575       U     df )       O           df,
576         I           myTime, myIter, myThid )
577  #endif  #endif
578         ELSE         ELSE
579          PRINT*,'invalid tracer indentity: ', tracerIdentity          PRINT*,'invalid tracer indentity: ', tracerIdentity
# Line 420  C *note* should update KPP_TRANSPORT_T t Line 581  C *note* should update KPP_TRANSPORT_T t
581         ENDIF         ENDIF
582         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
583          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
584           fVerT(i,j,kUp) = fVerT(i,j,kUp) + df(i,j)*maskUp(i,j)           fVerT(i,j,kUp) = fVerT(i,j,kUp)
585         &                  + df(i,j)*maskUp(i,j)*rhoFacF(k)
586          ENDDO          ENDDO
587         ENDDO         ENDDO
588        ENDIF        ENDIF
589  #endif  #endif
590    
591  C--   Divergence of fluxes  C--   Divergence of fluxes
592    C     Anelastic: scale vertical fluxes by rhoFac and leave Horizontal fluxes unchanged
593        DO j=1-Oly,sNy+Oly-1        DO j=1-Oly,sNy+Oly-1
594         DO i=1-Olx,sNx+Olx-1         DO i=1-Olx,sNx+Olx-1
595          gTracer(i,j,k,bi,bj)=gTracer(i,j,k,bi,bj)          gTracer(i,j,k,bi,bj)=gTracer(i,j,k,bi,bj)
596       &   -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &   -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
597       &    *recip_rA(i,j,bi,bj)       &   *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)*recip_rhoFacC(k)
598       &    *(       &   *( (fZon(i+1,j)-fZon(i,j))
599       &    +( fZon(i+1,j)-fZon(i,j) )       &     +(fMer(i,j+1)-fMer(i,j))
600       &    +( fMer(i,j+1)-fMer(i,j) )       &     +(fVerT(i,j,kDown)-fVerT(i,j,kUp))*rkSign
601       &    +( fVerT(i,j,kUp)-fVerT(i,j,kDown) )*rkFac       &     -localT(i,j)*( (uTrans(i+1,j)-uTrans(i,j))
602         &                   +(vTrans(i,j+1)-vTrans(i,j))
603         &                   +(rTransKp1(i,j)-rTrans(i,j))*rAdvFac
604         &                  )*advFac
605       &    )       &    )
606         ENDDO         ENDDO
607        ENDDO        ENDDO
608    
609  #ifdef NONLIN_FRSURF  #ifdef ALLOW_DEBUG
610  C-- account for 3.D divergence of the flow in rStar coordinate:        IF ( debugLevel .GE. debLevB
611        IF (calcAdvection .AND. select_rStar.GT.0) THEN       &   .AND. tracerIdentity.EQ.GAD_TEMPERATURE
612         DO j=1-Oly,sNy+Oly-1       &   .AND. k.EQ.2 .AND. myIter.EQ.1+nIter0
613          DO i=1-Olx,sNx+Olx-1       &   .AND. nPx.EQ.1 .AND. nPy.EQ.1
614           gTracer(i,j,k,bi,bj) = gTracer(i,j,k,bi,bj)       &   .AND. useCubedSphereExchange ) THEN
615       &     - (rStarExpC(i,j,bi,bj) - 1. _d 0)/deltaTfreesurf          CALL DEBUG_CS_CORNER_UV( ' fZon,fMer from GAD_CALC_RHS',
616       &       *tracer(i,j,k,bi,bj)*maskC(i,j,k,bi,bj)       &             fZon,fMer, k, standardMessageUnit,bi,bj,myThid )
         ENDDO  
        ENDDO  
       ENDIF  
       IF (calcAdvection .AND. select_rStar.LT.0) THEN  
        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  
617        ENDIF        ENDIF
618  #endif /* NONLIN_FRSURF */  #endif /* ALLOW_DEBUG */
         
619    
620        RETURN        RETURN
621        END        END
Legend:
Removed from v.1.19  
changed lines
  Added in v.1.47
  ViewVC Help
Powered by ViewVC 1.1.22