/[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.3 by adcroft, Thu Aug 30 00:40:37 2001 UTC revision 1.17 by jmc, Sun Jan 26 21:08:36 2003 UTC
# Line 3  C $Name$ Line 3  C $Name$
3    
4  #include "GAD_OPTIONS.h"  #include "GAD_OPTIONS.h"
5    
6    CBOP
7    C !ROUTINE: GAD_CALC_RHS
8    
9    C !INTERFACE: ==========================================================
10        SUBROUTINE GAD_CALC_RHS(        SUBROUTINE GAD_CALC_RHS(
11       I           bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,       I           bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,
12       I           xA,yA,uTrans,vTrans,rTrans,maskUp,       I           xA,yA,uTrans,vTrans,rTrans,maskUp,
13       I           diffKh, diffK4, KappaRT, Tracer,       I           diffKh, diffK4, KappaRT, Tracer,
14       I           tracerIdentity, advectionScheme,       I           tracerIdentity, advectionScheme, calcAdvection,
15       U           fVerT, gTracer,       U           fVerT, gTracer,
16       I           myThid )       I           myThid )
 C     /==========================================================\  
 C     | SUBROUTINE GAD_CALC_RHS                                  |  
 C     |==========================================================|  
 C     \==========================================================/  
       IMPLICIT NONE  
17    
18  C     == GLobal variables ==  C !DESCRIPTION:
19    C Calculates the tendancy of a tracer due to advection and diffusion.
20    C It calculates the fluxes in each direction indepentently and then
21    C sets the tendancy to the divergence of these fluxes. The advective
22    C fluxes are only calculated here when using the linear advection schemes
23    C otherwise only the diffusive and parameterized fluxes are calculated.
24    C
25    C Contributions to the flux are calculated and added:
26    C \begin{equation*}
27    C {\bf F} = {\bf F}_{adv} + {\bf F}_{diff} +{\bf F}_{GM} + {\bf F}_{KPP}
28    C \end{equation*}
29    C
30    C The tendancy is the divergence of the fluxes:
31    C \begin{equation*}
32    C G_\theta = G_\theta + \nabla \cdot {\bf F}
33    C \end{equation*}
34    C
35    C The tendancy is assumed to contain data on entry.
36    
37    C !USES: ===============================================================
38          IMPLICIT NONE
39  #include "SIZE.h"  #include "SIZE.h"
40  #include "EEPARAMS.h"  #include "EEPARAMS.h"
41  #include "PARAMS.h"  #include "PARAMS.h"
42  #include "GRID.h"  #include "GRID.h"
43  #include "DYNVARS.h"  #include "DYNVARS.h"
44    #include "SURFACE.h"
45  #include "GAD.h"  #include "GAD.h"
46    
47  C     == Routine arguments ==  #ifdef ALLOW_AUTODIFF_TAMC
48        INTEGER k,kUp,kDown,kM1  #include "tamc.h"
49    #include "tamc_keys.h"
50    #endif /* ALLOW_AUTODIFF_TAMC */
51    
52    C !INPUT PARAMETERS: ===================================================
53    C  bi,bj                :: tile indices
54    C  iMin,iMax,jMin,jMax  :: loop range for called routines
55    C  kup                  :: index into 2 1/2D array, toggles between 1 and 2
56    C  kdown                :: index into 2 1/2D array, toggles between 2 and 1
57    C  kp1                  :: =k+1 for k<Nr, =Nr for k=Nr
58    C  xA,yA                :: areas of X and Y face of tracer cells
59    C  uTrans,vTrans,rTrans :: 2-D arrays of volume transports at U,V and W points
60    C  maskUp               :: 2-D array for mask at W points
61    C  diffKh               :: horizontal diffusion coefficient
62    C  diffK4               :: bi-harmonic diffusion coefficient
63    C  KappaRT              :: 3-D array for vertical diffusion coefficient
64    C  Tracer               :: tracer field
65    C  tracerIdentity       :: identifier for the tracer (required only for KPP)
66    C  advectionScheme      :: advection scheme to use
67    C  calcAdvection        :: =False if Advec terms computed with multiDim scheme
68    C  myThid               :: thread number
69        INTEGER bi,bj,iMin,iMax,jMin,jMax        INTEGER bi,bj,iMin,iMax,jMin,jMax
70          INTEGER k,kUp,kDown,kM1
71        _RS xA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS xA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
72        _RS yA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS yA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
73        _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
# Line 38  C     == Routine arguments == Line 79  C     == Routine arguments ==
79        _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)
80        INTEGER tracerIdentity        INTEGER tracerIdentity
81        INTEGER advectionScheme        INTEGER advectionScheme
82        _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        LOGICAL calcAdvection
       _RL gTracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)  
83        INTEGER myThid        INTEGER myThid
84    
85  C     == Local variables ==  C !OUTPUT PARAMETERS: ==================================================
86  C     I, J, K - Loop counters  C  gTracer              :: tendancy array
87    C  fVerT                :: 2 1/2D arrays for vertical advective flux
88          _RL gTracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
89          _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
90    
91    C !LOCAL VARIABLES: ====================================================
92    C  i,j                  :: loop indices
93    C  df4                  :: used for storing del^2 T for bi-harmonic term
94    C  fZon                 :: zonal flux
95    C  fmer                 :: meridional flux
96    C  af                   :: advective flux
97    C  df                   :: diffusive flux
98    C  localT               :: local copy of tracer field
99        INTEGER i,j        INTEGER i,j
       LOGICAL TOP_LAYER  
       _RL afFacT, dfFacT  
100        _RL df4   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL df4   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
101        _RL fZon  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fZon  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
102        _RL fMer  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fMer  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
103        _RL af    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL af    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
104        _RL df    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL df    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
105        _RL localT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL localT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
106    CEOP
107    
108  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
109  C--   only the kUp part of fverT is set in this subroutine  C--   only the kUp part of fverT is set in this subroutine
110  C--   the kDown is still required  C--   the kDown is still required
111        fVerT(1,1,kDown) = fVerT(1,1,kDown)        fVerT(1,1,kDown) = fVerT(1,1,kDown)
112  #endif  #endif
113    
114        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
115         DO i=1-OLx,sNx+OLx         DO i=1-OLx,sNx+OLx
116          fZon(i,j)      = 0.0          fZon(i,j)      = 0. _d 0
117          fMer(i,j)      = 0.0          fMer(i,j)      = 0. _d 0
118          fVerT(i,j,kUp) = 0.0          fVerT(i,j,kUp) = 0. _d 0
119            df(i,j)        = 0. _d 0
120            df4(i,j)       = 0. _d 0
121            localT(i,j)    = 0. _d 0
122         ENDDO         ENDDO
123        ENDDO        ENDDO
124    
       afFacT = 1. _d 0  
       dfFacT = 1. _d 0  
       TOP_LAYER = K .EQ. 1  
   
125  C--   Make local copy of tracer array  C--   Make local copy of tracer array
126        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
127         DO i=1-OLx,sNx+OLx         DO i=1-OLx,sNx+OLx
# Line 78  C--   Make local copy of tracer array Line 129  C--   Make local copy of tracer array
129         ENDDO         ENDDO
130        ENDDO        ENDDO
131    
132    C--   Unless we have already calculated the advection terms we initialize
133    C     the tendency to zero.
134          IF (calcAdvection) THEN
135           DO j=1-Oly,sNy+Oly
136            DO i=1-Olx,sNx+Olx
137             gTracer(i,j,k,bi,bj)=0. _d 0
138            ENDDO
139           ENDDO
140          ENDIF
141    
142  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
143        IF (diffK4 .NE. 0.) THEN        IF (diffK4 .NE. 0.) THEN
# Line 89  C--   Pre-calculate del^2 T if bi-harmon Line 149  C--   Pre-calculate del^2 T if bi-harmon
149  C--   Initialize net flux in X direction  C--   Initialize net flux in X direction
150        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
151         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
152          fZon(i,j) = 0.          fZon(i,j) = 0. _d 0
153         ENDDO         ENDDO
154        ENDDO        ENDDO
155    
156  C-    Advective flux in X  C-    Advective flux in X
157          IF (calcAdvection) THEN
158        IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN        IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
159         CALL GAD_C2_ADV_X(bi,bj,k,uTrans,localT,af,myThid)         CALL GAD_C2_ADV_X(bi,bj,k,uTrans,localT,af,myThid)
160        ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN        ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
# Line 103  C-    Advective flux in X Line 164  C-    Advective flux in X
164         CALL GAD_U3_ADV_X(bi,bj,k,uTrans,localT,af,myThid)         CALL GAD_U3_ADV_X(bi,bj,k,uTrans,localT,af,myThid)
165        ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN        ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
166         CALL GAD_C4_ADV_X(bi,bj,k,uTrans,localT,af,myThid)         CALL GAD_C4_ADV_X(bi,bj,k,uTrans,localT,af,myThid)
167          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
168           CALL GAD_DST3_ADV_X(
169         &       bi,bj,k,deltaTtracer,uTrans,uVel,localT,af,myThid)
170          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
171           CALL GAD_DST3FL_ADV_X(
172         &       bi,bj,k,deltaTtracer,uTrans,uVel,localT,af,myThid)
173        ELSE        ELSE
174         STOP 'GAD_CALC_RHS: Bad advectionScheme (X)'         STOP 'GAD_CALC_RHS: Bad advectionScheme (X)'
175        ENDIF        ENDIF
176        DO j=jMin,jMax        DO j=1-Oly,sNy+Oly
177         DO i=iMin,iMax         DO i=1-Olx,sNx+Olx
178          fZon(i,j) = fZon(i,j) + af(i,j)          fZon(i,j) = fZon(i,j) + af(i,j)
179         ENDDO         ENDDO
180        ENDDO        ENDDO
181          ENDIF
182    
183  C-    Diffusive flux in X  C-    Diffusive flux in X
184        IF (diffKh.NE.0.) THEN        IF (diffKh.NE.0.) THEN
185         CALL GAD_DIFF_X(bi,bj,k,xA,diffKh,localT,df,myThid)         CALL GAD_DIFF_X(bi,bj,k,xA,diffKh,localT,df,myThid)
186        ELSE        ELSE
187         DO j=jMin,jMax         DO j=1-Oly,sNy+Oly
188          DO i=iMin,iMax          DO i=1-Olx,sNx+Olx
189           df(i,j) = 0.           df(i,j) = 0. _d 0
190          ENDDO          ENDDO
191         ENDDO         ENDDO
192        ENDIF        ENDIF
# Line 129  C-    GM/Redi flux in X Line 197  C-    GM/Redi flux in X
197  C *note* should update GMREDI_XTRANSPORT to use localT and set df  *aja*  C *note* should update GMREDI_XTRANSPORT to use localT and set df  *aja*
198          CALL GMREDI_XTRANSPORT(          CALL GMREDI_XTRANSPORT(
199       I     iMin,iMax,jMin,jMax,bi,bj,K,       I     iMin,iMax,jMin,jMax,bi,bj,K,
200       I     xA,Tracer,       I     xA,Tracer,tracerIdentity,
201       U     df,       U     df,
202       I     myThid)       I     myThid)
203        ENDIF        ENDIF
204  #endif  #endif
205        DO j=jMin,jMax        DO j=1-Oly,sNy+Oly
206         DO i=iMin,iMax         DO i=1-Olx,sNx+Olx
207          fZon(i,j) = fZon(i,j) + df(i,j)          fZon(i,j) = fZon(i,j) + df(i,j)
208         ENDDO         ENDDO
209        ENDDO        ENDDO
# Line 143  C *note* should update GMREDI_XTRANSPORT Line 211  C *note* should update GMREDI_XTRANSPORT
211  C-    Bi-harmonic duffusive flux in X  C-    Bi-harmonic duffusive flux in X
212        IF (diffK4 .NE. 0.) THEN        IF (diffK4 .NE. 0.) THEN
213         CALL GAD_BIHARM_X(bi,bj,k,xA,df4,diffK4,df,myThid)         CALL GAD_BIHARM_X(bi,bj,k,xA,df4,diffK4,df,myThid)
214         DO j=jMin,jMax         DO j=1-Oly,sNy+Oly
215          DO i=iMin,iMax          DO i=1-Olx,sNx+Olx
216           fZon(i,j) = fZon(i,j) + df(i,j)           fZon(i,j) = fZon(i,j) + df(i,j)
217          ENDDO          ENDDO
218         ENDDO         ENDDO
# Line 153  C-    Bi-harmonic duffusive flux in X Line 221  C-    Bi-harmonic duffusive flux in X
221  C--   Initialize net flux in Y direction  C--   Initialize net flux in Y direction
222        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
223         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
224          fMer(i,j) = 0.          fMer(i,j) = 0. _d 0
225         ENDDO         ENDDO
226        ENDDO        ENDDO
227    
228  C-    Advective flux in Y  C-    Advective flux in Y
229          IF (calcAdvection) THEN
230        IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN        IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
231         CALL GAD_C2_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)         CALL GAD_C2_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)
232        ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN        ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
# Line 167  C-    Advective flux in Y Line 236  C-    Advective flux in Y
236         CALL GAD_U3_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)         CALL GAD_U3_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)
237        ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN        ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
238         CALL GAD_C4_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)         CALL GAD_C4_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)
239          ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
240           CALL GAD_DST3_ADV_Y(
241         &       bi,bj,k,deltaTtracer,vTrans,vVel,localT,af,myThid)
242          ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
243           CALL GAD_DST3FL_ADV_Y(
244         &       bi,bj,k,deltaTtracer,vTrans,vVel,localT,af,myThid)
245        ELSE        ELSE
246         STOP 'GAD_CALC_RHS: Bad advectionScheme (Y)'         STOP 'GAD_CALC_RHS: Bad advectionScheme (Y)'
247        ENDIF        ENDIF
# Line 175  C-    Advective flux in Y Line 250  C-    Advective flux in Y
250          fMer(i,j) = fMer(i,j) + af(i,j)          fMer(i,j) = fMer(i,j) + af(i,j)
251         ENDDO         ENDDO
252        ENDDO        ENDDO
253          ENDIF
254    
255  C-    Diffusive flux in Y  C-    Diffusive flux in Y
256        IF (diffKh.NE.0.) THEN        IF (diffKh.NE.0.) THEN
# Line 182  C-    Diffusive flux in Y Line 258  C-    Diffusive flux in Y
258        ELSE        ELSE
259         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
260          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
261           df(i,j) = 0.           df(i,j) = 0. _d 0
262          ENDDO          ENDDO
263         ENDDO         ENDDO
264        ENDIF        ENDIF
# Line 190  C-    Diffusive flux in Y Line 266  C-    Diffusive flux in Y
266  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
267  C-    GM/Redi flux in Y  C-    GM/Redi flux in Y
268        IF (useGMRedi) THEN        IF (useGMRedi) THEN
        CALL GMREDI_YTRANSPORT(  
269  C *note* should update GMREDI_YTRANSPORT to use localT and set df  *aja*  C *note* should update GMREDI_YTRANSPORT to use localT and set df  *aja*
270           CALL GMREDI_YTRANSPORT(
271       I     iMin,iMax,jMin,jMax,bi,bj,K,       I     iMin,iMax,jMin,jMax,bi,bj,K,
272       I     yA,Tracer,       I     yA,Tracer,tracerIdentity,
273       U     df,       U     df,
274       I     myThid)       I     myThid)
275        ENDIF        ENDIF
# Line 214  C-    Bi-harmonic flux in Y Line 290  C-    Bi-harmonic flux in Y
290         ENDDO         ENDDO
291        ENDIF        ENDIF
292    
293  C--   Initialize net flux in R  #ifdef NONLIN_FRSURF
294        DO j=jMin,jMax  C--   Compute vertical flux fVerT(kDown) at interface k+1 (between k & k+1):
295         DO i=iMin,iMax        IF ( calcAdvection .AND. K.EQ.Nr .AND.
296          fVerT(i,j,kUp) = 0.       &     useRealFreshWaterFlux .AND.
297         ENDDO       &     buoyancyRelation .EQ. 'OCEANICP' ) THEN  
298        ENDDO         DO j=1-Oly,sNy+Oly
299            DO i=1-Olx,sNx+Olx
300             fVerT(i,j,kDown) = convertEmP2rUnit*PmEpR(i,j,bi,bj)
301         &     *rA(i,j,bi,bj)*maskC(i,j,k,bi,bj)*Tracer(i,j,k,bi,bj)
302            ENDDO
303           ENDDO
304          ENDIF
305    #endif /* NONLIN_FRSURF */
306    
307    C--   Compute vertical flux fVerT(kUp) at interface k (between k-1 & k):
308  C-    Advective flux in R  C-    Advective flux in R
309          IF (calcAdvection) THEN
310  C     Note: wVel needs to be masked  C     Note: wVel needs to be masked
311        IF (K.GE.2) THEN        IF (K.GE.2) THEN
312  C-    Compute vertical advective flux in the interior:  C-    Compute vertical advective flux in the interior:
# Line 234  C-    Compute vertical advective flux in Line 319  C-    Compute vertical advective flux in
319          CALL GAD_U3_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)          CALL GAD_U3_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)
320         ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN         ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
321          CALL GAD_C4_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)          CALL GAD_C4_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)
322  c       CALL GAD_C2_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)         ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
323            CALL GAD_DST3_ADV_R(
324         &       bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)
325           ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
326            CALL GAD_DST3FL_ADV_R(
327         &       bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)
328         ELSE         ELSE
329          STOP 'GAD_CALC_RHS: Bad advectionScheme (R)'          STOP 'GAD_CALC_RHS: Bad advectionScheme (R)'
330         ENDIF         ENDIF
# Line 255  C-    Surface "correction" term at k=1 : Line 345  C-    Surface "correction" term at k=1 :
345         ENDDO         ENDDO
346        ENDIF        ENDIF
347  C-    add the advective flux to fVerT  C-    add the advective flux to fVerT
348        DO j=jMin,jMax        DO j=1-Oly,sNy+Oly
349         DO i=iMin,iMax         DO i=1-Olx,sNx+Olx
350          fVerT(i,j,kUp) = fVerT(i,j,kUp) + afFacT*af(i,j)          fVerT(i,j,kUp) = fVerT(i,j,kUp) + af(i,j)
351         ENDDO         ENDDO
352        ENDDO        ENDDO
353          ENDIF
354    
355  C-    Diffusive flux in R  C-    Diffusive flux in R
356  C     Note: For K=1 then KM1=1 and this gives a dT/dr = 0 upper  C     Note: For K=1 then KM1=1 and this gives a dT/dr = 0 upper
357  C           boundary condition.  C           boundary condition.
358        IF (implicitDiffusion) THEN        IF (implicitDiffusion) THEN
359         DO j=jMin,jMax         DO j=1-Oly,sNy+Oly
360          DO i=iMin,iMax          DO i=1-Olx,sNx+Olx
361           df(i,j) = 0.           df(i,j) = 0. _d 0
362          ENDDO          ENDDO
363         ENDDO         ENDDO
364        ELSE        ELSE
365         CALL GAD_DIFF_R(bi,bj,k,KappaRT,tracer,df,myThid)         CALL GAD_DIFF_R(bi,bj,k,KappaRT,tracer,df,myThid)
366        ENDIF        ENDIF
 c     DO j=jMin,jMax  
 c      DO i=iMin,iMax  
 c       fVerT(i,j,kUp) = fVerT(i,j,kUp) + dfFacT*df(i,j)*maskUp(i,j)  
 c      ENDDO  
 c     ENDDO  
367    
368  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
369  C-    GM/Redi flux in R  C-    GM/Redi flux in R
# Line 285  C-    GM/Redi flux in R Line 371  C-    GM/Redi flux in R
371  C *note* should update GMREDI_RTRANSPORT to set df  *aja*  C *note* should update GMREDI_RTRANSPORT to set df  *aja*
372         CALL GMREDI_RTRANSPORT(         CALL GMREDI_RTRANSPORT(
373       I     iMin,iMax,jMin,jMax,bi,bj,K,       I     iMin,iMax,jMin,jMax,bi,bj,K,
374       I     maskUp,Tracer,       I     Tracer,tracerIdentity,
375       U     df,       U     df,
376       I     myThid)       I     myThid)
 c      DO j=jMin,jMax  
 c       DO i=iMin,iMax  
 c        fVerT(i,j,kUp) = fVerT(i,j,kUp) + dfFacT*df(i,j)*maskUp(i,j)  
 c       ENDDO  
 c      ENDDO  
377        ENDIF        ENDIF
378  #endif  #endif
379    
380        DO j=jMin,jMax        DO j=1-Oly,sNy+Oly
381         DO i=iMin,iMax         DO i=1-Olx,sNx+Olx
382          fVerT(i,j,kUp) = fVerT(i,j,kUp) + dfFacT*df(i,j)*maskUp(i,j)          fVerT(i,j,kUp) = fVerT(i,j,kUp) + df(i,j)*maskUp(i,j)
383         ENDDO         ENDDO
384        ENDDO        ENDDO
385    
386  #ifdef ALLOW_KPP  #ifdef ALLOW_KPP
387  C-    Add non local KPP transport term (ghat) to diffusive T flux.  C-    Add non local KPP transport term (ghat) to diffusive T flux.
388        IF (useKPP) THEN        IF (useKPP) THEN
389         DO j=jMin,jMax         DO j=1-Oly,sNy+Oly
390          DO i=iMin,iMax          DO i=1-Olx,sNx+Olx
391           df(i,j) = 0.           df(i,j) = 0. _d 0
392          ENDDO          ENDDO
393         ENDDO         ENDDO
394         IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN         IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
# Line 324  C *note* should update KPP_TRANSPORT_T t Line 405  C *note* should update KPP_TRANSPORT_T t
405         ELSE         ELSE
406          STOP 'GAD_CALC_RHS: Ooops'          STOP 'GAD_CALC_RHS: Ooops'
407         ENDIF         ENDIF
408         DO j=jMin,jMax         DO j=1-Oly,sNy+Oly
409          DO i=iMin,iMax          DO i=1-Olx,sNx+Olx
410           fVerT(i,j,kUp) = fVerT(i,j,kUp) + dfFacT*df(i,j)*maskUp(i,j)           fVerT(i,j,kUp) = fVerT(i,j,kUp) + df(i,j)*maskUp(i,j)
411          ENDDO          ENDDO
412         ENDDO         ENDDO
413        ENDIF        ENDIF
414  #endif  #endif
415    
416  C--   Divergence of fluxes  C--   Divergence of fluxes
417        DO j=jMin,jMax        DO j=1-Oly,sNy+Oly-1
418         DO i=iMin,iMax         DO i=1-Olx,sNx+Olx-1
419          gTracer(i,j,k,bi,bj)=          gTracer(i,j,k,bi,bj)=gTracer(i,j,k,bi,bj)
420       &   -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &   -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
421       &    *recip_rA(i,j,bi,bj)       &    *recip_rA(i,j,bi,bj)
422       &    *(       &    *(
# Line 346  C--   Divergence of fluxes Line 427  C--   Divergence of fluxes
427         ENDDO         ENDDO
428        ENDDO        ENDDO
429    
430    #ifdef NONLIN_FRSURF
431    C-- account for 3.D divergence of the flow in rStar coordinate:
432          IF (calcAdvection .AND. select_rStar.GT.0) THEN
433           DO j=1-Oly,sNy+Oly-1
434            DO i=1-Olx,sNx+Olx-1
435             gTracer(i,j,k,bi,bj) = gTracer(i,j,k,bi,bj)
436         &     - (rStarExpC(i,j,bi,bj) - 1. _d 0)/deltaTfreesurf
437         &       *tracer(i,j,k,bi,bj)*maskC(i,j,k,bi,bj)
438            ENDDO
439           ENDDO
440          ENDIF
441          IF (calcAdvection .AND. select_rStar.LT.0) THEN
442           DO j=1-Oly,sNy+Oly-1
443            DO i=1-Olx,sNx+Olx-1
444             gTracer(i,j,k,bi,bj) = gTracer(i,j,k,bi,bj)
445         &     - rStarDhCDt(i,j,bi,bj)
446         &       *tracer(i,j,k,bi,bj)*maskC(i,j,k,bi,bj)
447            ENDDO
448           ENDDO
449          ENDIF
450    #endif /* NONLIN_FRSURF */
451          
452    
453        RETURN        RETURN
454        END        END
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.17
  ViewVC Help
Powered by ViewVC 1.1.22