/[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.8 by jmc, Thu Dec 16 22:28:43 2004 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          INTEGER     kk
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 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)
      &     b5d(i,j,k) = -deltaTtracer  
137       &                  *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
# 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)
      &     d5d(i,j,k) = -deltaTtracer  
147       &                 *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
# Line 196  C--   Residual transp = Bolus transp + E Line 204  C--   Residual transp = Bolus transp + E
204            DO i=iMin,iMax            DO i=iMin,iMax
205  c          localTijk(i,j,k)     = gTracer(i,j,k,bi,bj)  c          localTijk(i,j,k)     = gTracer(i,j,k,bi,bj)
206             gTracer(i,j,k,bi,bj) = gTracer(i,j,k,bi,bj)             gTracer(i,j,k,bi,bj) = gTracer(i,j,k,bi,bj)
207       &      + deltaTtracer*recip_rA(i,j,bi,bj)       &      + dTtracerLev(1)*recip_rA(i,j,bi,bj)
208       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
209       &       *tracer(i,j,k,bi,bj)*(rTrans(i,j)-rTransKp1(i,j))*rkFac       &       *tracer(i,j,k,bi,bj)*(rTrans(i,j)-rTransKp1(i,j))*rkFac
210            ENDDO            ENDDO
211          ENDDO          ENDDO
212    
213          IF (K.GE.2) THEN  #ifdef ALLOW_AIM
214    C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr
215            IF ( K.GE.2 .AND.
216         &     (.NOT.useAIM .OR.tracerIdentity.NE.GAD_SALINITY .OR.K.LT.Nr)
217         &              ) THEN
218    #else
219            IF ( K.GE.2 ) THEN
220    #endif
221    
222           IF ( advectionScheme.EQ.ENUM_CENTERED_2ND ) THEN           IF ( advectionScheme.EQ.ENUM_CENTERED_2ND ) THEN
223            diagonalNumber = 3            diagonalNumber = 3
224            CALL GAD_C2_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,            CALL GAD_C2_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,
225       I                        deltaTtracer, rTrans,       I                        dTtracerLev, rTrans,
226       U                        b5d, c5d, d5d,       U                        b5d, c5d, d5d,
227       I                        myThid)       I                        myThid)
228           ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN           ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
229            diagonalNumber = 3            diagonalNumber = 3
230            CALL GAD_FLUXLIMIT_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,            CALL GAD_FLUXLIMIT_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,
231       I                        deltaTtracer, rTrans, localTijk,       I                        dTtracerLev, rTrans, localTijk,
232       U                        b5d, c5d, d5d,       U                        b5d, c5d, d5d,
233       I                        myThid)       I                        myThid)
234           ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD .OR.           ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD .OR.
235       &           advectionScheme.EQ.ENUM_CENTERED_4TH) THEN       &           advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
236            diagonalNumber = 5            diagonalNumber = 5
237            CALL GAD_U3C4_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,            CALL GAD_U3C4_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,
238       I                        advectionScheme, deltaTtracer, rTrans,       I                        advectionScheme, dTtracerLev, rTrans,
239       U                        a5d, b5d, c5d, d5d, e5d,       U                        a5d, b5d, c5d, d5d, e5d,
240       I                        myThid)       I                        myThid)
241           ELSE           ELSE
# Line 259  C--   Solve penta-diagonal system : Line 274  C--   Solve penta-diagonal system :
274          STOP 'GAD_IMPLICIT_R: no solver available'          STOP 'GAD_IMPLICIT_R: no solver available'
275        ENDIF        ENDIF
276    
277    #ifdef ALLOW_DIAGNOSTICS
278    C--   Set diagnostic suffix for the current tracer
279          IF ( useDiagnostics .AND. implicitDiffusion ) THEN
280            diagSufx = GAD_DIAG_SUFX( tracerIdentity, myThid )
281            diagName = 'DFrI'//diagSufx
282            IF ( DIAGNOSTICS_IS_ON(diagName,myThid) ) THEN
283             DO k= 1,Nr
284              IF ( k.EQ.1 ) THEN
285    C-  Note: Needs to call DIAGNOSTICS_FILL at level k=1 even if array == 0
286    C         otherwise counter is not incremented !!
287                DO j=1-OLy,sNy+OLy
288                 DO i=1-OLx,sNx+OLx
289                   df(i,j) = 0. _d 0
290                 ENDDO
291                ENDDO
292              ELSE
293                DO j=1,sNy
294                 DO i=1,sNx
295                   df(i,j) =
296         &              rA(i,j,bi,bj)
297         &            * KappaRX(i,j,k)*recip_drC(k)
298         &            * (gTracer(i,j,k,bi,bj) - gTracer(i,j,k-1,bi,bj))
299                 ENDDO
300                ENDDO
301              ENDIF
302              kk = -k
303              CALL DIAGNOSTICS_FILL(df,diagName, kk,1, 2,bi,bj, myThid)
304             ENDDO
305            ENDIF
306          ENDIF
307    #endif /* ALLOW_DIAGNOSTICS */
308    
309    C--   end if Nr > 1
310          ENDIF
311    
312        RETURN        RETURN
313        END        END
Legend:
Removed from v.1.4  
changed lines
  Added in v.1.8
  ViewVC Help
Powered by ViewVC 1.1.22