/[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.5 by jmc, Fri Nov 12 17:31:16 2004 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
79          CHARACTER*8 diagName
80          CHARACTER*4 GAD_DIAG_SUFX, diagSufx
81          EXTERNAL    GAD_DIAG_SUFX
82          LOGICAL     DIAGNOSTICS_IS_ON
83          EXTERNAL    DIAGNOSTICS_IS_ON
84          _RL df (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
85    #endif
86  CEOP  CEOP
87    
88        IF (Nr.LE.1) RETURN  C--   no need to solve anything with only 1 level:
89          IF (Nr.GT.1) THEN
90    
91  C--   Initialise  C--   Initialise
92        iMin = 1        iMin = 1
# Line 95  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 123  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            IF (maskC(i,j,k-1,bi,bj).EQ.1.)             b5d(i,j,k) = -deltaTLev(k)*maskC(i,j,k-1,bi,bj)
139       &     b5d(i,j,k) = -deltaTtracer       &                  *_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 134  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            IF (maskC(i,j,k+1,bi,bj).EQ.1.)             d5d(i,j,k) = -deltaTLev(k)*maskC(i,j,k+1,bi,bj)
149       &     d5d(i,j,k) = -deltaTtracer       &                 *_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 175  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       &      + deltaTtracer*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)*(rTrans(i,j)-rTransKp1(i,j))*rkFac       &       *tracer(i,j,k,bi,bj)*(rTransKp1(i,j)-rTrans(i,j))*rkSign
214            ENDDO            ENDDO
215          ENDDO          ENDDO
216    
# Line 214  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                        deltaTtracer, rTrans,       I                        deltaTLev, rTrans,
230         U                        b5d, c5d, d5d,
231         I                        myThid )
232             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,       U                        b5d, c5d, d5d,
238       I                        myThid)       I                        myThid )
239           ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN           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                        deltaTtracer, 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, deltaTtracer, rTrans,       I                        advectionScheme, deltaTLev, rTrans,
251       U                        a5d, b5d, c5d, d5d, e5d,       U                        a5d, b5d, c5d, d5d, e5d,
252       I                        myThid)       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,
258         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 266  C--   Solve penta-diagonal system : Line 292  C--   Solve penta-diagonal system :
292          STOP 'GAD_IMPLICIT_R: no solver available'          STOP 'GAD_IMPLICIT_R: no solver available'
293        ENDIF        ENDIF
294    
295    #ifdef ALLOW_DIAGNOSTICS
296    C--   Set diagnostic suffix for the current tracer
297          IF ( useDiagnostics .AND. implicitDiffusion ) THEN
298            diagSufx = GAD_DIAG_SUFX( tracerIdentity, myThid )
299            diagName = 'DFrI'//diagSufx
300            IF ( DIAGNOSTICS_IS_ON(diagName,myThid) ) THEN
301             DO k= 1,Nr
302              IF ( k.EQ.1 ) THEN
303    C-  Note: Needs to call DIAGNOSTICS_FILL at level k=1 even if array == 0
304    C         otherwise counter is not incremented !!
305                DO j=1-OLy,sNy+OLy
306                 DO i=1-OLx,sNx+OLx
307                   df(i,j) = 0. _d 0
308                 ENDDO
309                ENDDO
310              ELSE
311                DO j=1,sNy
312                 DO i=1,sNx
313                   df(i,j) =
314         &              rA(i,j,bi,bj)
315         &            * KappaRX(i,j,k)*recip_drC(k)
316         &            * (gTracer(i,j,k,bi,bj) - gTracer(i,j,k-1,bi,bj))
317                 ENDDO
318                ENDDO
319              ENDIF
320              CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
321             ENDDO
322            ENDIF
323          ENDIF
324    #endif /* ALLOW_DIAGNOSTICS */
325    
326    C--   end if Nr > 1
327          ENDIF
328    
329        RETURN        RETURN
330        END        END
Legend:
Removed from v.1.5  
changed lines
  Added in v.1.14
  ViewVC Help
Powered by ViewVC 1.1.22