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

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

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

revision 1.10 by jmc, Wed Jun 22 00:27:47 2005 UTC revision 1.14 by jahn, Fri Jun 26 23:10:09 2009 UTC
# Line 6  C $Name$ Line 6  C $Name$
6  CBOP  CBOP
7  C     !ROUTINE: GAD_IMPLICIT_R  C     !ROUTINE: GAD_IMPLICIT_R
8  C     !INTERFACE:  C     !INTERFACE:
9        SUBROUTINE GAD_IMPLICIT_R(        SUBROUTINE GAD_IMPLICIT_R(
10       I      implicitAdvection, advectionScheme, tracerIdentity,       I      implicitAdvection, advectionScheme, tracerIdentity,
11       I      kappaRX, wVel, tracer,       I      deltaTLev,
12         I      kappaRX, wVel, tracer,
13       U      gTracer,       U      gTracer,
14       I      bi, bj, myTime, myIter, myThid )       I      bi, bj, myTime, myIter, myThid )
15  C     !DESCRIPTION:  C     !DESCRIPTION:
# Line 39  C myThid            :: thread number Line 40  C myThid            :: thread number
40        LOGICAL implicitAdvection        LOGICAL implicitAdvection
41        INTEGER advectionScheme        INTEGER advectionScheme
42        INTEGER tracerIdentity        INTEGER tracerIdentity
43          _RL     deltaTLev(Nr)
44        _RL kappaRX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL kappaRX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
45        _RL wVel   (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)        _RL wVel   (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
46        _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)
# Line 69  C errCode   :: > 0 if singular matrix Line 71  C errCode   :: > 0 if singular matrix
71        _RL c5d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL c5d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
72        _RL d5d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL d5d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
73        _RL e5d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL e5d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
74        _RL rTrans(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL wFld     (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
75          _RL rTrans   (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
76        _RL rTransKp1(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL rTransKp1(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
77        _RL localTijk(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL localTijk(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
78  #ifdef ALLOW_DIAGNOSTICS  #ifdef ALLOW_DIAGNOSTICS
# Line 104  C--   Initialise Line 107  C--   Initialise
107        diagonalNumber = 1        diagonalNumber = 1
108    
109  C--   Non-Linear Advection scheme: keep a local copy of tracer field  C--   Non-Linear Advection scheme: keep a local copy of tracer field
110        IF ( advectionScheme.EQ.ENUM_FLUX_LIMIT .OR.        IF ( advectionScheme.EQ.ENUM_FLUX_LIMIT .OR.
111       &     advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN       &     advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
112          IF ( multiDimAdvection ) THEN          IF ( multiDimAdvection ) THEN
113           DO k=1,Nr           DO k=1,Nr
# Line 132  C-     1rst lower diagonal : Line 135  C-     1rst lower diagonal :
135         DO k=2,Nr         DO k=2,Nr
136          DO j=jMin,jMax          DO j=jMin,jMax
137           DO i=iMin,iMax           DO i=iMin,iMax
138             b5d(i,j,k) = -dTtracerLev(k)*maskC(i,j,k-1,bi,bj)             b5d(i,j,k) = -deltaTLev(k)*maskC(i,j,k-1,bi,bj)
139       &                  *recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &                  *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
140       &                  *kappaRX(i,j, k )*recip_drC( k )       &                  *kappaRX(i,j, k )*recip_drC( k )
141           ENDDO           ENDDO
142          ENDDO          ENDDO
# Line 142  C-     1rst upper diagonal : Line 145  C-     1rst upper diagonal :
145         DO k=1,Nr-1         DO k=1,Nr-1
146          DO j=jMin,jMax          DO j=jMin,jMax
147           DO i=iMin,iMax           DO i=iMin,iMax
148             d5d(i,j,k) = -dTtracerLev(k)*maskC(i,j,k+1,bi,bj)             d5d(i,j,k) = -deltaTLev(k)*maskC(i,j,k+1,bi,bj)
149       &                 *recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &                 *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
150       &                 *KappaRX(i,j,k+1)*recip_drC(k+1)       &                 *KappaRX(i,j,k+1)*recip_drC(k+1)
151           ENDDO           ENDDO
152          ENDDO          ENDDO
# Line 182  C--    Compute transport Line 185  C--    Compute transport
185          IF (k.EQ.1) THEN          IF (k.EQ.1) THEN
186           DO j=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
187            DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
188              rTrans(i,j) = 0.              wFld(i,j)   = 0. _d 0
189                rTrans(i,j) = 0. _d 0
190            ENDDO            ENDDO
191           ENDDO           ENDDO
192          ELSE          ELSE
193           DO j=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
194            DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
195              rTrans(i,j) = wVel(i,j,k,bi,bj)*rA(i,j,bi,bj)              wFld(i,j)   = wVel(i,j,k,bi,bj)
196       &                   *maskC(i,j,k-1,bi,bj)              rTrans(i,j) = wFld(i,j)*rA(i,j,bi,bj)*maskC(i,j,k-1,bi,bj)
197            ENDDO            ENDDO
198           ENDDO           ENDDO
199  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
200  C--   Residual transp = Bolus transp + Eulerian transp  C--   Residual transp = Bolus transp + Eulerian transp
201            IF (useGMRedi)           IF (useGMRedi)
202       &      CALL GMREDI_CALC_WFLOW(       &     CALL GMREDI_CALC_WFLOW(
203       &                    rTrans, bi, bj, k, myThid)       U                 wFld, rTrans,
204         I                 k, bi, bj, myThid )
205  #endif /* ALLOW_GMREDI */  #endif /* ALLOW_GMREDI */
206          ENDIF          ENDIF
207          DO j=jMin,jMax          DO j=jMin,jMax
208            DO i=iMin,iMax            DO i=iMin,iMax
209  c          localTijk(i,j,k)     = gTracer(i,j,k,bi,bj)  c          localTijk(i,j,k)     = gTracer(i,j,k,bi,bj)
210             gTracer(i,j,k,bi,bj) = gTracer(i,j,k,bi,bj)             gTracer(i,j,k,bi,bj) = gTracer(i,j,k,bi,bj)
211       &      + dTtracerLev(1)*recip_rA(i,j,bi,bj)       &      + deltaTLev(1)*recip_rA(i,j,bi,bj)
212       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
213       &       *tracer(i,j,k,bi,bj)*(rTransKp1(i,j)-rTrans(i,j))*rkSign       &       *tracer(i,j,k,bi,bj)*(rTransKp1(i,j)-rTrans(i,j))*rkSign
214            ENDDO            ENDDO
# Line 221  C- a hack to prevent Water-Vapor vert.tr Line 226  C- a hack to prevent Water-Vapor vert.tr
226           IF ( advectionScheme.EQ.ENUM_CENTERED_2ND ) THEN           IF ( advectionScheme.EQ.ENUM_CENTERED_2ND ) THEN
227            diagonalNumber = 3            diagonalNumber = 3
228            CALL GAD_C2_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,            CALL GAD_C2_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,
229       I                        dTtracerLev, rTrans,       I                        deltaTLev, rTrans,
230       U                        b5d, c5d, d5d,       U                        b5d, c5d, d5d,
231       I                        myThid)       I                        myThid )
232           ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN           ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST
233         &       .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
234              diagonalNumber = 3
235              CALL GAD_DST2U1_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,
236         I                        advectionScheme, deltaTLev, rTrans,
237         U                        b5d, c5d, d5d,
238         I                        myThid )
239             ELSEIF ( advectionScheme.EQ.ENUM_FLUX_LIMIT ) THEN
240            diagonalNumber = 3            diagonalNumber = 3
241            CALL GAD_FLUXLIMIT_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,            CALL GAD_FLUXLIMIT_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,
242       I                        dTtracerLev, rTrans, localTijk,       I                        deltaTLev, rTrans, localTijk,
243       U                        b5d, c5d, d5d,       U                        b5d, c5d, d5d,
244       I                        myThid)       I                        myThid )
245           ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD .OR.           ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_3RD
246       &           advectionScheme.EQ.ENUM_CENTERED_4TH) THEN       &       .OR. advectionScheme.EQ.ENUM_CENTERED_4TH
247         &       .OR. advectionScheme.EQ.ENUM_DST3 ) THEN
248            diagonalNumber = 5            diagonalNumber = 5
249            CALL GAD_U3C4_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,            CALL GAD_U3C4_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,
250       I                        advectionScheme, dTtracerLev, rTrans,       I                        advectionScheme, deltaTLev, rTrans,
251         U                        a5d, b5d, c5d, d5d, e5d,
252         I                        myThid )
253             ELSEIF ( advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
254              diagonalNumber = 5
255              CALL GAD_DST3FL_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,
256         I                        deltaTLev, rTrans, localTijk,
257       U                        a5d, b5d, c5d, d5d, e5d,       U                        a5d, b5d, c5d, d5d, e5d,
258       I                        myThid)       I                        myThid )
259           ELSE           ELSE
260            STOP 'GAD_IMPLICIT_R: Adv.Scheme in Impl form not yet coded'            STOP 'GAD_IMPLICIT_R: Adv.Scheme in Impl form not yet coded'
261           ENDIF           ENDIF
# Line 291  C         otherwise counter is not incre Line 310  C         otherwise counter is not incre
310            ELSE            ELSE
311              DO j=1,sNy              DO j=1,sNy
312               DO i=1,sNx               DO i=1,sNx
313                 df(i,j) =                 df(i,j) =
314       &              rA(i,j,bi,bj)       &              rA(i,j,bi,bj)
315       &            * KappaRX(i,j,k)*recip_drC(k)       &            * KappaRX(i,j,k)*recip_drC(k)
316       &            * (gTracer(i,j,k,bi,bj) - gTracer(i,j,k-1,bi,bj))       &            * (gTracer(i,j,k,bi,bj) - gTracer(i,j,k-1,bi,bj))
Legend:
Removed from v.1.10  
changed lines
  Added in v.1.14
  ViewVC Help
Powered by ViewVC 1.1.22