/[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.10 by jmc, Wed Jun 22 00:27:47 2005 UTC
# Line 72  C errCode   :: > 0 if singular matrix Line 72  C errCode   :: > 0 if singular matrix
72        _RL rTrans(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL rTrans(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
73        _RL rTransKp1(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL rTransKp1(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
74        _RL localTijk(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL localTijk(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
75    #ifdef ALLOW_DIAGNOSTICS
76          CHARACTER*8 diagName
77          CHARACTER*4 GAD_DIAG_SUFX, diagSufx
78          EXTERNAL    GAD_DIAG_SUFX
79          LOGICAL     DIAGNOSTICS_IS_ON
80          EXTERNAL    DIAGNOSTICS_IS_ON
81          _RL df (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
82    #endif
83  CEOP  CEOP
84    
85        IF (Nr.LE.1) RETURN  C--   no need to solve anything with only 1 level:
86          IF (Nr.GT.1) THEN
87    
88  C--   Initialise  C--   Initialise
89        iMin = 1        iMin = 1
# Line 123  C-     1rst lower diagonal : Line 132  C-     1rst lower diagonal :
132         DO k=2,Nr         DO k=2,Nr
133          DO j=jMin,jMax          DO j=jMin,jMax
134           DO i=iMin,iMax           DO i=iMin,iMax
135            IF (maskC(i,j,k-1,bi,bj).EQ.1.)             b5d(i,j,k) = -dTtracerLev(k)*maskC(i,j,k-1,bi,bj)
      &     b5d(i,j,k) = -deltaTtracer  
136       &                  *recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &                  *recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
137       &                  *kappaRX(i,j, k )*recip_drC( k )       &                  *kappaRX(i,j, k )*recip_drC( k )
138           ENDDO           ENDDO
# Line 134  C-     1rst upper diagonal : Line 142  C-     1rst upper diagonal :
142         DO k=1,Nr-1         DO k=1,Nr-1
143          DO j=jMin,jMax          DO j=jMin,jMax
144           DO i=iMin,iMax           DO i=iMin,iMax
145            IF (maskC(i,j,k+1,bi,bj).EQ.1.)             d5d(i,j,k) = -dTtracerLev(k)*maskC(i,j,k+1,bi,bj)
      &     d5d(i,j,k) = -deltaTtracer  
146       &                 *recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &                 *recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
147       &                 *KappaRX(i,j,k+1)*recip_drC(k+1)       &                 *KappaRX(i,j,k+1)*recip_drC(k+1)
148           ENDDO           ENDDO
# Line 196  C--   Residual transp = Bolus transp + E Line 203  C--   Residual transp = Bolus transp + E
203            DO i=iMin,iMax            DO i=iMin,iMax
204  c          localTijk(i,j,k)     = gTracer(i,j,k,bi,bj)  c          localTijk(i,j,k)     = gTracer(i,j,k,bi,bj)
205             gTracer(i,j,k,bi,bj) = gTracer(i,j,k,bi,bj)             gTracer(i,j,k,bi,bj) = gTracer(i,j,k,bi,bj)
206       &      + deltaTtracer*recip_rA(i,j,bi,bj)       &      + dTtracerLev(1)*recip_rA(i,j,bi,bj)
207       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
208       &       *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
209            ENDDO            ENDDO
210          ENDDO          ENDDO
211    
# Line 214  C- a hack to prevent Water-Vapor vert.tr Line 221  C- a hack to prevent Water-Vapor vert.tr
221           IF ( advectionScheme.EQ.ENUM_CENTERED_2ND ) THEN           IF ( advectionScheme.EQ.ENUM_CENTERED_2ND ) THEN
222            diagonalNumber = 3            diagonalNumber = 3
223            CALL GAD_C2_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,            CALL GAD_C2_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,
224       I                        deltaTtracer, rTrans,       I                        dTtracerLev, rTrans,
225       U                        b5d, c5d, d5d,       U                        b5d, c5d, d5d,
226       I                        myThid)       I                        myThid)
227           ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN           ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
228            diagonalNumber = 3            diagonalNumber = 3
229            CALL GAD_FLUXLIMIT_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,            CALL GAD_FLUXLIMIT_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,
230       I                        deltaTtracer, rTrans, localTijk,       I                        dTtracerLev, rTrans, localTijk,
231       U                        b5d, c5d, d5d,       U                        b5d, c5d, d5d,
232       I                        myThid)       I                        myThid)
233           ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD .OR.           ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD .OR.
234       &           advectionScheme.EQ.ENUM_CENTERED_4TH) THEN       &           advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
235            diagonalNumber = 5            diagonalNumber = 5
236            CALL GAD_U3C4_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,            CALL GAD_U3C4_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,
237       I                        advectionScheme, deltaTtracer, rTrans,       I                        advectionScheme, dTtracerLev, rTrans,
238       U                        a5d, b5d, c5d, d5d, e5d,       U                        a5d, b5d, c5d, d5d, e5d,
239       I                        myThid)       I                        myThid)
240           ELSE           ELSE
# Line 266  C--   Solve penta-diagonal system : Line 273  C--   Solve penta-diagonal system :
273          STOP 'GAD_IMPLICIT_R: no solver available'          STOP 'GAD_IMPLICIT_R: no solver available'
274        ENDIF        ENDIF
275    
276    #ifdef ALLOW_DIAGNOSTICS
277    C--   Set diagnostic suffix for the current tracer
278          IF ( useDiagnostics .AND. implicitDiffusion ) THEN
279            diagSufx = GAD_DIAG_SUFX( tracerIdentity, myThid )
280            diagName = 'DFrI'//diagSufx
281            IF ( DIAGNOSTICS_IS_ON(diagName,myThid) ) THEN
282             DO k= 1,Nr
283              IF ( k.EQ.1 ) THEN
284    C-  Note: Needs to call DIAGNOSTICS_FILL at level k=1 even if array == 0
285    C         otherwise counter is not incremented !!
286                DO j=1-OLy,sNy+OLy
287                 DO i=1-OLx,sNx+OLx
288                   df(i,j) = 0. _d 0
289                 ENDDO
290                ENDDO
291              ELSE
292                DO j=1,sNy
293                 DO i=1,sNx
294                   df(i,j) =
295         &              rA(i,j,bi,bj)
296         &            * KappaRX(i,j,k)*recip_drC(k)
297         &            * (gTracer(i,j,k,bi,bj) - gTracer(i,j,k-1,bi,bj))
298                 ENDDO
299                ENDDO
300              ENDIF
301              CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
302             ENDDO
303            ENDIF
304          ENDIF
305    #endif /* ALLOW_DIAGNOSTICS */
306    
307    C--   end if Nr > 1
308          ENDIF
309    
310        RETURN        RETURN
311        END        END
Legend:
Removed from v.1.5  
changed lines
  Added in v.1.10
  ViewVC Help
Powered by ViewVC 1.1.22