/[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.5 by adcroft, Tue Sep 4 17:00:48 2001 UTC revision 1.25 by jmc, Fri Jun 25 18:19:20 2004 UTC
# Line 3  C $Name$ Line 3  C $Name$
3    
4  #include "GAD_OPTIONS.h"  #include "GAD_OPTIONS.h"
5    
6    CBOP
7    C !ROUTINE: GAD_CALC_RHS
8    
9    C !INTERFACE: ==========================================================
10        SUBROUTINE GAD_CALC_RHS(        SUBROUTINE GAD_CALC_RHS(
11       I           bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,       I           bi,bj,iMin,iMax,jMin,jMax,k,kM1,kUp,kDown,
12       I           xA,yA,uTrans,vTrans,rTrans,maskUp,       I           xA,yA,uTrans,vTrans,rTrans,rTransKp1,maskUp,
13         I           uVel, vVel, wVel,
14       I           diffKh, diffK4, KappaRT, Tracer,       I           diffKh, diffK4, KappaRT, Tracer,
15       I           tracerIdentity, advectionScheme,       I           tracerIdentity, advectionScheme,
16         I           calcAdvection, implicitAdvection,
17       U           fVerT, gTracer,       U           fVerT, gTracer,
18       I           myThid )       I           myThid )
 C     /==========================================================\  
 C     | SUBROUTINE GAD_CALC_RHS                                  |  
 C     |==========================================================|  
 C     \==========================================================/  
       IMPLICIT NONE  
19    
20  C     == GLobal variables ==  C !DESCRIPTION:
21    C Calculates the tendancy of a tracer due to advection and diffusion.
22    C It calculates the fluxes in each direction indepentently and then
23    C sets the tendancy to the divergence of these fluxes. The advective
24    C fluxes are only calculated here when using the linear advection schemes
25    C otherwise only the diffusive and parameterized fluxes are calculated.
26    C
27    C Contributions to the flux are calculated and added:
28    C \begin{equation*}
29    C {\bf F} = {\bf F}_{adv} + {\bf F}_{diff} +{\bf F}_{GM} + {\bf F}_{KPP}
30    C \end{equation*}
31    C
32    C The tendancy is the divergence of the fluxes:
33    C \begin{equation*}
34    C G_\theta = G_\theta + \nabla \cdot {\bf F}
35    C \end{equation*}
36    C
37    C The tendancy is assumed to contain data on entry.
38    
39    C !USES: ===============================================================
40          IMPLICIT NONE
41  #include "SIZE.h"  #include "SIZE.h"
42  #include "EEPARAMS.h"  #include "EEPARAMS.h"
43  #include "PARAMS.h"  #include "PARAMS.h"
44  #include "GRID.h"  #include "GRID.h"
45  #include "DYNVARS.h"  #include "SURFACE.h"
46  #include "GAD.h"  #include "GAD.h"
47    
48  C     == Routine arguments ==  #ifdef ALLOW_AUTODIFF_TAMC
49        INTEGER k,kUp,kDown,kM1  #include "tamc.h"
50    #include "tamc_keys.h"
51    #endif /* ALLOW_AUTODIFF_TAMC */
52    
53    C !INPUT PARAMETERS: ===================================================
54    C bi,bj            :: tile indices
55    C iMin,iMax        :: loop range for called routines
56    C jMin,jMax        :: loop range for called routines
57    C kup              :: index into 2 1/2D array, toggles between 1|2
58    C kdown            :: index into 2 1/2D array, toggles between 2|1
59    C kp1              :: =k+1 for k<Nr, =Nr for k=Nr
60    C xA,yA            :: areas of X and Y face of tracer cells
61    C uTrans,vTrans    :: 2-D arrays of volume transports at U,V points
62    C rTrans           :: 2-D arrays of volume transports at W points
63    C rTransKp1        :: 2-D array of volume trans at W pts, interf k+1
64    C maskUp           :: 2-D array for mask at W points
65    C uVel,vVel,wVel   :: 3 components of the velcity field (3-D array)
66    C diffKh           :: horizontal diffusion coefficient
67    C diffK4           :: bi-harmonic diffusion coefficient
68    C KappaRT          :: 3-D array for vertical diffusion coefficient
69    C Tracer           :: tracer field
70    C tracerIdentity   :: tracer identifier (required for KPP,GM)
71    C advectionScheme  :: advection scheme to use
72    C calcAdvection    :: =False if Advec computed with multiDim scheme
73    C implicitAdvection:: =True if vertical Advec computed implicitly
74    C myThid           :: thread number
75        INTEGER bi,bj,iMin,iMax,jMin,jMax        INTEGER bi,bj,iMin,iMax,jMin,jMax
76          INTEGER k,kUp,kDown,kM1
77        _RS xA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS xA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
78        _RS yA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS yA    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
79        _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
80        _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
81        _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
82          _RL rTransKp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
83        _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
84          _RL uVel  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
85          _RL vVel  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
86          _RL wVel  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
87        _RL diffKh, diffK4        _RL diffKh, diffK4
88        _RL KappaRT(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL KappaRT(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
89        _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)
90        INTEGER tracerIdentity        INTEGER tracerIdentity
91        INTEGER advectionScheme        INTEGER advectionScheme
92        _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)        LOGICAL calcAdvection
93        _RL gTracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        LOGICAL implicitAdvection
94        INTEGER myThid        INTEGER myThid
95    
96  C     == Local variables ==  C !OUTPUT PARAMETERS: ==================================================
97  C     I, J, K - Loop counters  C gTracer          :: tendancy array
98    C fVerT            :: 2 1/2D arrays for vertical advective flux
99          _RL gTracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
100          _RL fVerT (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
101    
102    C !LOCAL VARIABLES: ====================================================
103    C i,j              :: loop indices
104    C df4              :: used for storing del^2 T for bi-harmonic term
105    C fZon             :: zonal flux
106    C fmer             :: meridional flux
107    C af               :: advective flux
108    C df               :: diffusive flux
109    C localT           :: local copy of tracer field
110        INTEGER i,j        INTEGER i,j
       LOGICAL TOP_LAYER  
       _RL afFacT, dfFacT  
111        _RL df4   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL df4   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
112        _RL fZon  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fZon  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
113        _RL fMer  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL fMer  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
114        _RL af    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL af    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
115        _RL df    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL df    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
116        _RL localT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL localT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
117          _RL advFac, rAdvFac
118    CEOP
119    
120  #ifdef ALLOW_AUTODIFF_TAMC  #ifdef ALLOW_AUTODIFF_TAMC
121  C--   only the kUp part of fverT is set in this subroutine  C--   only the kUp part of fverT is set in this subroutine
122  C--   the kDown is still required  C--   the kDown is still required
123        fVerT(1,1,kDown) = fVerT(1,1,kDown)        fVerT(1,1,kDown) = fVerT(1,1,kDown)
124  #endif  #endif
125    
126          advFac  = 0. _d 0
127          IF (calcAdvection) advFac = 1. _d 0
128          rAdvFac = rkFac*advFac
129          IF (implicitAdvection) rAdvFac = 0. _d 0
130    
131        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
132         DO i=1-OLx,sNx+OLx         DO i=1-OLx,sNx+OLx
133          fZon(i,j)      = 0.0          fZon(i,j)      = 0. _d 0
134          fMer(i,j)      = 0.0          fMer(i,j)      = 0. _d 0
135          fVerT(i,j,kUp) = 0.0          fVerT(i,j,kUp) = 0. _d 0
136            df(i,j)        = 0. _d 0
137            df4(i,j)       = 0. _d 0
138         ENDDO         ENDDO
139        ENDDO        ENDDO
140    
       afFacT = 1. _d 0  
       dfFacT = 1. _d 0  
       TOP_LAYER = K .EQ. 1  
   
141  C--   Make local copy of tracer array  C--   Make local copy of tracer array
142        DO j=1-OLy,sNy+OLy        DO j=1-OLy,sNy+OLy
143         DO i=1-OLx,sNx+OLx         DO i=1-OLx,sNx+OLx
# Line 78  C--   Make local copy of tracer array Line 145  C--   Make local copy of tracer array
145         ENDDO         ENDDO
146        ENDDO        ENDDO
147    
148    C--   Unless we have already calculated the advection terms we initialize
149    C     the tendency to zero.
150    C     <== now done earlier at the beginning of thermodynamics.
151    c     IF (calcAdvection) THEN
152    c      DO j=1-Oly,sNy+Oly
153    c       DO i=1-Olx,sNx+Olx
154    c        gTracer(i,j,k,bi,bj)=0. _d 0
155    c       ENDDO
156    c      ENDDO
157    c     ENDIF
158    
159  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
160        IF (diffK4 .NE. 0.) THEN        IF (diffK4 .NE. 0.) THEN
# Line 89  C--   Pre-calculate del^2 T if bi-harmon Line 166  C--   Pre-calculate del^2 T if bi-harmon
166  C--   Initialize net flux in X direction  C--   Initialize net flux in X direction
167        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
168         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
169          fZon(i,j) = 0.          fZon(i,j) = 0. _d 0
170         ENDDO         ENDDO
171        ENDDO        ENDDO
172    
173  C-    Advective flux in X  C-    Advective flux in X
174          IF (calcAdvection) THEN
175        IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN        IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
176         CALL GAD_C2_ADV_X(bi,bj,k,uTrans,localT,af,myThid)         CALL GAD_C2_ADV_X(bi,bj,k,uTrans,localT,af,myThid)
177        ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN        ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
# Line 117  C-    Advective flux in X Line 195  C-    Advective flux in X
195          fZon(i,j) = fZon(i,j) + af(i,j)          fZon(i,j) = fZon(i,j) + af(i,j)
196         ENDDO         ENDDO
197        ENDDO        ENDDO
198          ENDIF
199    
200  C-    Diffusive flux in X  C-    Diffusive flux in X
201        IF (diffKh.NE.0.) THEN        IF (diffKh.NE.0.) THEN
# Line 124  C-    Diffusive flux in X Line 203  C-    Diffusive flux in X
203        ELSE        ELSE
204         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
205          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
206           df(i,j) = 0.           df(i,j) = 0. _d 0
207          ENDDO          ENDDO
208         ENDDO         ENDDO
209        ENDIF        ENDIF
# Line 135  C-    GM/Redi flux in X Line 214  C-    GM/Redi flux in X
214  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*
215          CALL GMREDI_XTRANSPORT(          CALL GMREDI_XTRANSPORT(
216       I     iMin,iMax,jMin,jMax,bi,bj,K,       I     iMin,iMax,jMin,jMax,bi,bj,K,
217       I     xA,Tracer,       I     xA,Tracer,tracerIdentity,
218       U     df,       U     df,
219       I     myThid)       I     myThid)
220        ENDIF        ENDIF
# Line 159  C-    Bi-harmonic duffusive flux in X Line 238  C-    Bi-harmonic duffusive flux in X
238  C--   Initialize net flux in Y direction  C--   Initialize net flux in Y direction
239        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
240         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
241          fMer(i,j) = 0.          fMer(i,j) = 0. _d 0
242         ENDDO         ENDDO
243        ENDDO        ENDDO
244    
245  C-    Advective flux in Y  C-    Advective flux in Y
246          IF (calcAdvection) THEN
247        IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN        IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
248         CALL GAD_C2_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)         CALL GAD_C2_ADV_Y(bi,bj,k,vTrans,localT,af,myThid)
249        ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN        ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
# Line 187  C-    Advective flux in Y Line 267  C-    Advective flux in Y
267          fMer(i,j) = fMer(i,j) + af(i,j)          fMer(i,j) = fMer(i,j) + af(i,j)
268         ENDDO         ENDDO
269        ENDDO        ENDDO
270          ENDIF
271    
272  C-    Diffusive flux in Y  C-    Diffusive flux in Y
273        IF (diffKh.NE.0.) THEN        IF (diffKh.NE.0.) THEN
# Line 194  C-    Diffusive flux in Y Line 275  C-    Diffusive flux in Y
275        ELSE        ELSE
276         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
277          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
278           df(i,j) = 0.           df(i,j) = 0. _d 0
279          ENDDO          ENDDO
280         ENDDO         ENDDO
281        ENDIF        ENDIF
# Line 202  C-    Diffusive flux in Y Line 283  C-    Diffusive flux in Y
283  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
284  C-    GM/Redi flux in Y  C-    GM/Redi flux in Y
285        IF (useGMRedi) THEN        IF (useGMRedi) THEN
        CALL GMREDI_YTRANSPORT(  
286  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*
287           CALL GMREDI_YTRANSPORT(
288       I     iMin,iMax,jMin,jMax,bi,bj,K,       I     iMin,iMax,jMin,jMax,bi,bj,K,
289       I     yA,Tracer,       I     yA,Tracer,tracerIdentity,
290       U     df,       U     df,
291       I     myThid)       I     myThid)
292        ENDIF        ENDIF
# Line 226  C-    Bi-harmonic flux in Y Line 307  C-    Bi-harmonic flux in Y
307         ENDDO         ENDDO
308        ENDIF        ENDIF
309    
310  C--   Initialize net flux in R  C--   Compute vertical flux fVerT(kUp) at interface k (between k-1 & k):
       DO j=1-Oly,sNy+Oly  
        DO i=1-Olx,sNx+Olx  
         fVerT(i,j,kUp) = 0.  
        ENDDO  
       ENDDO  
   
311  C-    Advective flux in R  C-    Advective flux in R
312  C     Note: wVel needs to be masked  #ifdef ALLOW_AIM
313        IF (K.GE.2) THEN  C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr
314          IF (calcAdvection .AND. .NOT.implicitAdvection .AND. K.GE.2 .AND.
315         &     (.NOT.useAIM .OR.tracerIdentity.NE.GAD_SALINITY .OR.K.LT.Nr)
316         &   ) THEN
317    #else
318          IF (calcAdvection .AND. .NOT.implicitAdvection .AND. K.GE.2) THEN
319    #endif
320  C-    Compute vertical advective flux in the interior:  C-    Compute vertical advective flux in the interior:
321         IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN         IF (advectionScheme.EQ.ENUM_CENTERED_2ND) THEN
322          CALL GAD_C2_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)          CALL GAD_C2_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)
# Line 247  C-    Compute vertical advective flux in Line 328  C-    Compute vertical advective flux in
328         ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN         ELSEIF (advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
329          CALL GAD_C4_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)          CALL GAD_C4_ADV_R(bi,bj,k,rTrans,tracer,af,myThid)
330         ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN         ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
331  c       CALL GAD_DST3_ADV_R(          CALL GAD_DST3_ADV_R(
332  c    &       bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)       &       bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)
         STOP 'GAD_CALC_RHS: GAD_DST3_ADV_R not coded yet'  
333         ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN         ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
334  c       CALL GAD_DST3FL_ADV_R(          CALL GAD_DST3FL_ADV_R(
335  c    &       bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)       &       bi,bj,k,deltaTtracer,rTrans,wVel,tracer,af,myThid)
         STOP 'GAD_CALC_RHS: GAD_DST3FL_ADV_R not coded yet'  
336         ELSE         ELSE
337          STOP 'GAD_CALC_RHS: Bad advectionScheme (R)'          STOP 'GAD_CALC_RHS: Bad advectionScheme (R)'
338         ENDIF         ENDIF
339  C-    Surface "correction" term at k>1 :  C-     add the advective flux to fVerT
340         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
341          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
342           af(i,j) = af(i,j)           fVerT(i,j,kUp) = fVerT(i,j,kUp) + af(i,j)
      &           + (maskC(i,j,k,bi,bj)-maskC(i,j,k-1,bi,bj))*  
      &             rTrans(i,j)*Tracer(i,j,k,bi,bj)  
343          ENDDO          ENDDO
        ENDDO  
       ELSE  
 C-    Surface "correction" term at k=1 :  
        DO j=1-Oly,sNy+Oly  
         DO i=1-Olx,sNx+Olx  
          af(i,j) = rTrans(i,j)*Tracer(i,j,k,bi,bj)  
         ENDDO  
        ENDDO  
       ENDIF  
 C-    add the advective flux to fVerT  
       DO j=1-Oly,sNy+Oly  
        DO i=1-Olx,sNx+Olx  
         fVerT(i,j,kUp) = fVerT(i,j,kUp) + afFacT*af(i,j)  
344         ENDDO         ENDDO
345        ENDDO        ENDIF
346    
347  C-    Diffusive flux in R  C-    Diffusive flux in R
348  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
# Line 286  C           boundary condition. Line 350  C           boundary condition.
350        IF (implicitDiffusion) THEN        IF (implicitDiffusion) THEN
351         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
352          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
353           df(i,j) = 0.           df(i,j) = 0. _d 0
354          ENDDO          ENDDO
355         ENDDO         ENDDO
356        ELSE        ELSE
357         CALL GAD_DIFF_R(bi,bj,k,KappaRT,tracer,df,myThid)         CALL GAD_DIFF_R(bi,bj,k,KappaRT,tracer,df,myThid)
358        ENDIF        ENDIF
 c     DO j=1-Oly,sNy+Oly  
 c      DO i=1-Olx,sNx+Olx  
 c       fVerT(i,j,kUp) = fVerT(i,j,kUp) + dfFacT*df(i,j)*maskUp(i,j)  
 c      ENDDO  
 c     ENDDO  
359    
360  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
361  C-    GM/Redi flux in R  C-    GM/Redi flux in R
# Line 304  C-    GM/Redi flux in R Line 363  C-    GM/Redi flux in R
363  C *note* should update GMREDI_RTRANSPORT to set df  *aja*  C *note* should update GMREDI_RTRANSPORT to set df  *aja*
364         CALL GMREDI_RTRANSPORT(         CALL GMREDI_RTRANSPORT(
365       I     iMin,iMax,jMin,jMax,bi,bj,K,       I     iMin,iMax,jMin,jMax,bi,bj,K,
366       I     maskUp,Tracer,       I     Tracer,tracerIdentity,
367       U     df,       U     df,
368       I     myThid)       I     myThid)
 c      DO j=1-Oly,sNy+Oly  
 c       DO i=1-Olx,sNx+Olx  
 c        fVerT(i,j,kUp) = fVerT(i,j,kUp) + dfFacT*df(i,j)*maskUp(i,j)  
 c       ENDDO  
 c      ENDDO  
369        ENDIF        ENDIF
370  #endif  #endif
371    
372        DO j=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly
373         DO i=1-Olx,sNx+Olx         DO i=1-Olx,sNx+Olx
374          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)
375         ENDDO         ENDDO
376        ENDDO        ENDDO
377    
# Line 326  C-    Add non local KPP transport term ( Line 380  C-    Add non local KPP transport term (
380        IF (useKPP) THEN        IF (useKPP) THEN
381         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
382          DO i=1-Olx,sNx+Olx          DO i=1-Olx,sNx+Olx
383           df(i,j) = 0.           df(i,j) = 0. _d 0
384          ENDDO          ENDDO
385         ENDDO         ENDDO
386         IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN         IF (tracerIdentity.EQ.GAD_TEMPERATURE) THEN
# Line 340  C *note* should update KPP_TRANSPORT_T t Line 394  C *note* should update KPP_TRANSPORT_T t
394       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,       I     iMin,iMax,jMin,jMax,bi,bj,k,km1,
395       I     KappaRT,       I     KappaRT,
396       U     df )       U     df )
397    #ifdef ALLOW_PTRACERS
398           ELSEIF (tracerIdentity .GE. GAD_TR1) THEN
399            CALL KPP_TRANSPORT_PTR(
400         I     iMin,iMax,jMin,jMax,bi,bj,k,km1,
401         I     tracerIdentity-GAD_TR1+1,KappaRT,
402         U     df )
403    #endif
404         ELSE         ELSE
405            PRINT*,'invalid tracer indentity: ', tracerIdentity
406          STOP 'GAD_CALC_RHS: Ooops'          STOP 'GAD_CALC_RHS: Ooops'
407         ENDIF         ENDIF
408         DO j=1-Oly,sNy+Oly         DO j=1-Oly,sNy+Oly
409          DO i=1-Olx,sNx+Olx          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=1-Oly,sNy+Oly        DO j=1-Oly,sNy+Oly-1
418         DO i=1-Olx,sNx+Olx         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)*recip_rA(i,j,bi,bj)
421       &    *recip_rA(i,j,bi,bj)       &   *( (fZon(i+1,j)-fZon(i,j))
422       &    *(       &     +(fMer(i,j+1)-fMer(i,j))
423       &    +( fZon(i+1,j)-fZon(i,j) )       &     +(fVerT(i,j,kUp)-fVerT(i,j,kDown))*rkFac
424       &    +( fMer(i,j+1)-fMer(i,j) )       &     -localT(i,j)*( (uTrans(i+1,j)-uTrans(i,j))
425       &    +( fVerT(i,j,kUp)-fVerT(i,j,kDown) )*rkFac       &                   +(vTrans(i,j+1)-vTrans(i,j))
426         &                   +(rTrans(i,j)-rTransKp1(i,j))*rAdvFac
427         &                  )*advFac
428       &    )       &    )
429         ENDDO         ENDDO
430        ENDDO        ENDDO
Legend:
Removed from v.1.5  
changed lines
  Added in v.1.25
  ViewVC Help
Powered by ViewVC 1.1.22