/[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.1 by jmc, Sat Jan 3 00:43:01 2004 UTC revision 1.2 by jmc, Wed Jan 7 21:35:00 2004 UTC
# Line 53  C     == Local variables == Line 53  C     == Local variables ==
53        _RL d5d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL d5d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
54        _RL e5d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)        _RL e5d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
55        _RL rTrans(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL rTrans(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
56          _RL rTransKp1(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
57          _RL localTijk(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
58          _RL limitDf  (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
59        _RL rFlx        _RL rFlx
60  CEOP  CEOP
61    
# Line 76  C--   Initialise Line 79  C--   Initialise
79        ENDDO        ENDDO
80        diagonalNumber = 1        diagonalNumber = 1
81    
82    C--   Non-Linear Advection scheme: keep a local copy of tracer field
83          IF ( advectionScheme.EQ.ENUM_FLUX_LIMIT .OR.
84         &     advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
85            IF ( multiDimAdvection ) THEN
86             DO k=1,Nr
87              DO j=1-Oly,sNy+Oly
88               DO i=1-Olx,sNx+Olx
89                localTijk(i,j,k) = gTracer(i,j,k,bi,bj)
90               ENDDO
91              ENDDO
92             ENDDO
93            ELSE
94             DO k=1,Nr
95              DO j=1-Oly,sNy+Oly
96               DO i=1-Olx,sNx+Olx
97                localTijk(i,j,k) = tracer(i,j,k,bi,bj)
98               ENDDO
99              ENDDO
100             ENDDO
101            ENDIF
102          ENDIF
103    
104        IF (implicitDiffusion) THEN        IF (implicitDiffusion) THEN
105  C--   set the tri-diagonal matrix to solve the implicit diffusion problem  C--   set the tri-diagonal matrix to solve the implicit diffusion problem
106         diagonalNumber = 3         diagonalNumber = 3
# Line 115  C--   end if implicitDiffusion Line 140  C--   end if implicitDiffusion
140    
141        IF (implicitAdvection) THEN        IF (implicitAdvection) THEN
142    
143         IF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN         DO k=Nr,1,-1
          diagonalNumber = 3  
 C note: a) this should go into a separated gad_ S/R  
          DO k=2,Nr  
144    
145            DO j=1-Oly,sNy+Oly  C--    Compute transport
146             DO i=1-Olx,sNx+Olx          IF (k.EQ.Nr) THEN
147             DO j=1-Oly,sNy+Oly
148              DO i=1-Olx,sNx+Olx
149                rTransKp1(i,j) = 0.
150              ENDDO
151             ENDDO
152            ELSE
153             DO j=1-Oly,sNy+Oly
154              DO i=1-Olx,sNx+Olx
155                rTransKp1(i,j) = rTrans(i,j)
156              ENDDO
157             ENDDO
158            ENDIF
159    
160            IF (k.EQ.1) THEN
161             DO j=1-Oly,sNy+Oly
162              DO i=1-Olx,sNx+Olx
163                rTrans(i,j) = 0.
164              ENDDO
165             ENDDO
166            ELSE
167             DO j=1-Oly,sNy+Oly
168              DO i=1-Olx,sNx+Olx
169              rTrans(i,j) = wVel(i,j,k,bi,bj)*rA(i,j,bi,bj)              rTrans(i,j) = wVel(i,j,k,bi,bj)*rA(i,j,bi,bj)
170       &                   *maskC(i,j,k-1,bi,bj)       &                   *maskC(i,j,k-1,bi,bj)
            ENDDO  
171            ENDDO            ENDDO
172             ENDDO
173  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
174  C--   Residual transp = Bolus transp + Eulerian transp  C--   Residual transp = Bolus transp + Eulerian transp
175            IF (useGMRedi)            IF (useGMRedi)
176       &      CALL GMREDI_CALC_WFLOW(       &      CALL GMREDI_CALC_WFLOW(
177       &                    rTrans, bi, bj, k, myThid)       &                    rTrans, bi, bj, k, myThid)
178  #endif /* ALLOW_GMREDI */  #endif /* ALLOW_GMREDI */
179            ENDIF
180  C-       space Centered advection scheme, Flux form:          DO j=jMin,jMax
181            DO j=jMin,jMax            DO i=iMin,iMax
182             DO i=iMin,iMax  c          localTijk(i,j,k)     = gTracer(i,j,k,bi,bj)
183              rFlx = 0.5 _d 0 *deltaTtracer*rTrans(i,j)             gTracer(i,j,k,bi,bj) = gTracer(i,j,k,bi,bj)
184       &                      *recip_rA(i,j,bi,bj)*rkFac       &      + deltaTtracer*recip_rA(i,j,bi,bj)
185              b5d(i,j,k) = b5d(i,j,k)       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
186       &                 + rFlx*recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &       *tracer(i,j,k,bi,bj)*(rTrans(i,j)-rTransKp1(i,j))*rkFac
             c5d(i,j,k) = c5d(i,j,k)  
      &                 + rFlx*recip_hFacC(i,j,k,bi,bj)*recip_drF(k)  
             c5d(i,j,k-1) = c5d(i,j,k-1)  
      &                 - rFlx*recip_hFacC(i,j,k-1,bi,bj)*recip_drF(k-1)  
             d5d(i,j,k-1) = d5d(i,j,k-1)  
      &                 - rFlx*recip_hFacC(i,j,k-1,bi,bj)*recip_drF(k-1)  
            ENDDO  
187            ENDDO            ENDDO
188            ENDDO
189    
190  C--      end k loop          IF (K.GE.2) THEN
191           ENDDO  
192         ELSE           IF ( advectionScheme.EQ.ENUM_CENTERED_2ND ) THEN
193           STOP 'GAD_IMPLICIT_R: Adv.Scheme in Impl form not yet coded'            diagonalNumber = 3
194         ENDIF            CALL GAD_C2_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,
195         I                        deltaTtracer, rTrans,
196         U                        b5d, c5d, d5d,
197         I                        myThid)
198             ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
199              diagonalNumber = 3
200              CALL GAD_FLUXLIMIT_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,
201         I                        deltaTtracer, rTrans, localTijk,
202         U                        b5d, c5d, d5d,
203         I                        myThid)
204             ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD .OR.
205         &           advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
206              diagonalNumber = 5
207              CALL GAD_U3C4_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,
208         I                        advectionScheme, deltaTtracer, rTrans,
209         U                        a5d, b5d, c5d, d5d, e5d,
210         I                        myThid)
211             ELSE
212              STOP 'GAD_IMPLICIT_R: Adv.Scheme in Impl form not yet coded'
213             ENDIF
214    
215            ENDIF
216    
217    C--     end k loop
218           ENDDO
219    
220  C--   end if implicitAdvection  C--   end if implicitAdvection
221        ENDIF        ENDIF
# Line 168  C--   Solve tri-diagonal system : Line 230  C--   Solve tri-diagonal system :
230          IF (errCode.GE.1) THEN          IF (errCode.GE.1) THEN
231            STOP 'GAD_IMPLICIT_R: error when solving 3-Diag problem'            STOP 'GAD_IMPLICIT_R: error when solving 3-Diag problem'
232          ENDIF          ENDIF
233          ELSEIF ( diagonalNumber .EQ. 5 ) THEN
234    C--   Solve penta-diagonal system :
235            CALL SOLVE_PENTADIAGONAL( iMin,iMax, jMin,jMax,
236         I                            a5d, b5d, c5d, d5d, e5d,
237         U                            gTracer,
238         O                            errCode,
239         I                            bi, bj, myThid )
240            IF (errCode.GE.1) THEN
241              STOP 'GAD_IMPLICIT_R: error when solving 5-Diag problem'
242            ENDIF
243        ELSEIF ( diagonalNumber .NE. 1 ) THEN        ELSEIF ( diagonalNumber .NE. 1 ) THEN
244          STOP 'GAD_IMPLICIT_R: no solver available'          STOP 'GAD_IMPLICIT_R: no solver available'
245        ENDIF        ENDIF
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
  ViewVC Help
Powered by ViewVC 1.1.22