/[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.4 by edhill, Mon Mar 29 03:33:51 2004 UTC revision 1.13 by jmc, Sun Jun 18 23:27:44 2006 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      kappaRX, wVel, tracer,
12       U      gTracer,       U      gTracer,
13       I      bi, bj, myTime, myIter, myThid )       I      bi, bj, myTime, myIter, myThid )
14  C     !DESCRIPTION:  C     !DESCRIPTION:
# Line 69  C errCode   :: > 0 if singular matrix Line 69  C errCode   :: > 0 if singular matrix
69        _RL c5d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL c5d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
70        _RL d5d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL d5d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
71        _RL e5d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL e5d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
72        _RL rTrans(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL wFld     (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
73          _RL rTrans   (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
74        _RL rTransKp1(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL rTransKp1(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
75        _RL localTijk(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL localTijk(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
76    #ifdef ALLOW_DIAGNOSTICS
77          CHARACTER*8 diagName
78          CHARACTER*4 GAD_DIAG_SUFX, diagSufx
79          EXTERNAL    GAD_DIAG_SUFX
80          LOGICAL     DIAGNOSTICS_IS_ON
81          EXTERNAL    DIAGNOSTICS_IS_ON
82          _RL df (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
83    #endif
84  CEOP  CEOP
85    
86        IF (Nr.LE.1) RETURN  C--   no need to solve anything with only 1 level:
87          IF (Nr.GT.1) THEN
88    
89  C--   Initialise  C--   Initialise
90        iMin = 1        iMin = 1
# Line 95  C--   Initialise Line 105  C--   Initialise
105        diagonalNumber = 1        diagonalNumber = 1
106    
107  C--   Non-Linear Advection scheme: keep a local copy of tracer field  C--   Non-Linear Advection scheme: keep a local copy of tracer field
108        IF ( advectionScheme.EQ.ENUM_FLUX_LIMIT .OR.        IF ( advectionScheme.EQ.ENUM_FLUX_LIMIT .OR.
109       &     advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN       &     advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
110          IF ( multiDimAdvection ) THEN          IF ( multiDimAdvection ) THEN
111           DO k=1,Nr           DO k=1,Nr
# Line 123  C-     1rst lower diagonal : Line 133  C-     1rst lower diagonal :
133         DO k=2,Nr         DO k=2,Nr
134          DO j=jMin,jMax          DO j=jMin,jMax
135           DO i=iMin,iMax           DO i=iMin,iMax
136            IF (maskC(i,j,k-1,bi,bj).EQ.1.)             b5d(i,j,k) = -dTtracerLev(k)*maskC(i,j,k-1,bi,bj)
137       &     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)  
138       &                  *kappaRX(i,j, k )*recip_drC( k )       &                  *kappaRX(i,j, k )*recip_drC( k )
139           ENDDO           ENDDO
140          ENDDO          ENDDO
# Line 134  C-     1rst upper diagonal : Line 143  C-     1rst upper diagonal :
143         DO k=1,Nr-1         DO k=1,Nr-1
144          DO j=jMin,jMax          DO j=jMin,jMax
145           DO i=iMin,iMax           DO i=iMin,iMax
146            IF (maskC(i,j,k+1,bi,bj).EQ.1.)             d5d(i,j,k) = -dTtracerLev(k)*maskC(i,j,k+1,bi,bj)
147       &     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)  
148       &                 *KappaRX(i,j,k+1)*recip_drC(k+1)       &                 *KappaRX(i,j,k+1)*recip_drC(k+1)
149           ENDDO           ENDDO
150          ENDDO          ENDDO
# Line 175  C--    Compute transport Line 183  C--    Compute transport
183          IF (k.EQ.1) THEN          IF (k.EQ.1) THEN
184           DO j=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
185            DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
186              rTrans(i,j) = 0.              wFld(i,j)   = 0. _d 0
187                rTrans(i,j) = 0. _d 0
188            ENDDO            ENDDO
189           ENDDO           ENDDO
190          ELSE          ELSE
191           DO j=1-Oly,sNy+Oly           DO j=1-Oly,sNy+Oly
192            DO i=1-Olx,sNx+Olx            DO i=1-Olx,sNx+Olx
193              rTrans(i,j) = wVel(i,j,k,bi,bj)*rA(i,j,bi,bj)              wFld(i,j)   = wVel(i,j,k,bi,bj)
194       &                   *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)
195            ENDDO            ENDDO
196           ENDDO           ENDDO
197  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
198  C--   Residual transp = Bolus transp + Eulerian transp  C--   Residual transp = Bolus transp + Eulerian transp
199            IF (useGMRedi)           IF (useGMRedi)
200       &      CALL GMREDI_CALC_WFLOW(       &     CALL GMREDI_CALC_WFLOW(
201       &                    rTrans, bi, bj, k, myThid)       U                 wFld, rTrans,
202         I                 k, bi, bj, myThid )
203  #endif /* ALLOW_GMREDI */  #endif /* ALLOW_GMREDI */
204          ENDIF          ENDIF
205          DO j=jMin,jMax          DO j=jMin,jMax
206            DO i=iMin,iMax            DO i=iMin,iMax
207  c          localTijk(i,j,k)     = gTracer(i,j,k,bi,bj)  c          localTijk(i,j,k)     = gTracer(i,j,k,bi,bj)
208             gTracer(i,j,k,bi,bj) = gTracer(i,j,k,bi,bj)             gTracer(i,j,k,bi,bj) = gTracer(i,j,k,bi,bj)
209       &      + deltaTtracer*recip_rA(i,j,bi,bj)       &      + dTtracerLev(1)*recip_rA(i,j,bi,bj)
210       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
211       &       *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
212            ENDDO            ENDDO
213          ENDDO          ENDDO
214    
215          IF (K.GE.2) THEN  #ifdef ALLOW_AIM
216    C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr
217            IF ( K.GE.2 .AND.
218         &     (.NOT.useAIM .OR.tracerIdentity.NE.GAD_SALINITY .OR.K.LT.Nr)
219         &              ) THEN
220    #else
221            IF ( K.GE.2 ) THEN
222    #endif
223    
224           IF ( advectionScheme.EQ.ENUM_CENTERED_2ND ) THEN           IF ( advectionScheme.EQ.ENUM_CENTERED_2ND ) THEN
225            diagonalNumber = 3            diagonalNumber = 3
226            CALL GAD_C2_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,            CALL GAD_C2_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,
227       I                        deltaTtracer, rTrans,       I                        dTtracerLev, rTrans,
228       U                        b5d, c5d, d5d,       U                        b5d, c5d, d5d,
229       I                        myThid)       I                        myThid )
230           ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN           ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST
231         &       .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
232              diagonalNumber = 3
233              CALL GAD_DST2U1_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,
234         I                        advectionScheme, dTtracerLev, rTrans,
235         U                        b5d, c5d, d5d,
236         I                        myThid )
237             ELSEIF ( advectionScheme.EQ.ENUM_FLUX_LIMIT ) THEN
238            diagonalNumber = 3            diagonalNumber = 3
239            CALL GAD_FLUXLIMIT_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,            CALL GAD_FLUXLIMIT_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,
240       I                        deltaTtracer, rTrans, localTijk,       I                        dTtracerLev, rTrans, localTijk,
241       U                        b5d, c5d, d5d,       U                        b5d, c5d, d5d,
242       I                        myThid)       I                        myThid )
243           ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD .OR.           ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_3RD
244       &           advectionScheme.EQ.ENUM_CENTERED_4TH) THEN       &       .OR. advectionScheme.EQ.ENUM_CENTERED_4TH
245         &       .OR. advectionScheme.EQ.ENUM_DST3 ) THEN
246            diagonalNumber = 5            diagonalNumber = 5
247            CALL GAD_U3C4_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,            CALL GAD_U3C4_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,
248       I                        advectionScheme, deltaTtracer, rTrans,       I                        advectionScheme, dTtracerLev, rTrans,
249       U                        a5d, b5d, c5d, d5d, e5d,       U                        a5d, b5d, c5d, d5d, e5d,
250       I                        myThid)       I                        myThid )
251             ELSEIF ( advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
252              diagonalNumber = 5
253              CALL GAD_DST3FL_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,
254         I                        dTtracerLev, rTrans, localTijk,
255         U                        a5d, b5d, c5d, d5d, e5d,
256         I                        myThid )
257           ELSE           ELSE
258            STOP 'GAD_IMPLICIT_R: Adv.Scheme in Impl form not yet coded'            STOP 'GAD_IMPLICIT_R: Adv.Scheme in Impl form not yet coded'
259           ENDIF           ENDIF
# Line 259  C--   Solve penta-diagonal system : Line 290  C--   Solve penta-diagonal system :
290          STOP 'GAD_IMPLICIT_R: no solver available'          STOP 'GAD_IMPLICIT_R: no solver available'
291        ENDIF        ENDIF
292    
293    #ifdef ALLOW_DIAGNOSTICS
294    C--   Set diagnostic suffix for the current tracer
295          IF ( useDiagnostics .AND. implicitDiffusion ) THEN
296            diagSufx = GAD_DIAG_SUFX( tracerIdentity, myThid )
297            diagName = 'DFrI'//diagSufx
298            IF ( DIAGNOSTICS_IS_ON(diagName,myThid) ) THEN
299             DO k= 1,Nr
300              IF ( k.EQ.1 ) THEN
301    C-  Note: Needs to call DIAGNOSTICS_FILL at level k=1 even if array == 0
302    C         otherwise counter is not incremented !!
303                DO j=1-OLy,sNy+OLy
304                 DO i=1-OLx,sNx+OLx
305                   df(i,j) = 0. _d 0
306                 ENDDO
307                ENDDO
308              ELSE
309                DO j=1,sNy
310                 DO i=1,sNx
311                   df(i,j) =
312         &              rA(i,j,bi,bj)
313         &            * KappaRX(i,j,k)*recip_drC(k)
314         &            * (gTracer(i,j,k,bi,bj) - gTracer(i,j,k-1,bi,bj))
315                 ENDDO
316                ENDDO
317              ENDIF
318              CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
319             ENDDO
320            ENDIF
321          ENDIF
322    #endif /* ALLOW_DIAGNOSTICS */
323    
324    C--   end if Nr > 1
325          ENDIF
326    
327        RETURN        RETURN
328        END        END
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.13
  ViewVC Help
Powered by ViewVC 1.1.22