/[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.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: \bv  C     !DESCRIPTION:
16  C     *==========================================================*  C     Solve implicitly vertical advection and diffusion for one tracer.
 C     | S/R GAD_IMPLICIT_R  
 C     | o Solve implicitly vertical advection & diffusion  
 C     |   for one tracer  
 C     *==========================================================*  
 C     *==========================================================*  
 C     \ev  
17    
18  C     !USES:  C     !USES:
19        IMPLICIT NONE        IMPLICIT NONE
# Line 29  C     == Global data == Line 24  C     == Global data ==
24  #include "GRID.h"  #include "GRID.h"
25  #include "GAD.h"  #include "GAD.h"
26    
27  C     !INPUT/OUTPUT PARAMETERS:  C !INPUT/OUTPUT PARAMETERS:
28  C     == Routine Arguments ==  C == Routine Arguments ==
29    C implicitAdvection :: if True, treat vertical advection implicitly
30    C advectionScheme   :: advection scheme to use
31    C tracerIdentity    :: Identifier for the tracer
32    C kappaRX           :: 3-D array for vertical diffusion coefficient
33    C wVel              :: vertical component of the velcity field
34    C tracer            :: tracer field at current time step
35    C gTracer           :: future tracer field
36    C bi,bj             :: tile indices
37    C myTime            :: current time
38    C myIter            :: current iteration number
39    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 42  C     == Routine Arguments == Line 49  C     == Routine Arguments ==
49        _RL myTime        _RL myTime
50        INTEGER myIter, myThid        INTEGER myIter, myThid
51    
52  C     !LOCAL VARIABLES:  C !LOCAL VARIABLES:
53  C     == Local variables ==  C == Local variables ==
54    C iMin,iMax,jMin,jMax  :: computational domain
55    C i,j,k  :: loop indices
56    C a5d    :: 2nd  lower diagonal of the pentadiagonal matrix
57    C b5d    :: 1rst lower diagonal of the pentadiagonal matrix
58    C c5d    :: main diagonal       of the pentadiagonal matrix
59    C d5d    :: 1rst upper diagonal of the pentadiagonal matrix
60    C e5d    :: 2nd  upper diagonal of the pentadiagonal matrix
61    C rTrans    :: vertical volume transport at inteface k
62    C rTransKp1 :: vertical volume transport at inteface k+1
63    C localTijk :: local copy of tracer (for Non-Lin Adv.Scheme)
64    C diagonalNumber :: number of non-zero diagonals in the matrix
65    C errCode   :: > 0 if singular matrix
66        INTEGER iMin,iMax,jMin,jMax        INTEGER iMin,iMax,jMin,jMax
67        INTEGER i,j,k        INTEGER i,j,k
68        INTEGER diagonalNumber, errCode        INTEGER diagonalNumber, errCode
# Line 52  C     == Local variables == Line 71  C     == Local variables ==
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 rFlx        _RL rTrans   (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
76          _RL rTransKp1(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
77          _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 76  C--   Initialise Line 106  C--   Initialise
106        ENDDO        ENDDO
107        diagonalNumber = 1        diagonalNumber = 1
108    
109    C--   Non-Linear Advection scheme: keep a local copy of tracer field
110          IF ( advectionScheme.EQ.ENUM_FLUX_LIMIT .OR.
111         &     advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
112            IF ( multiDimAdvection ) THEN
113             DO k=1,Nr
114              DO j=1-Oly,sNy+Oly
115               DO i=1-Olx,sNx+Olx
116                localTijk(i,j,k) = gTracer(i,j,k,bi,bj)
117               ENDDO
118              ENDDO
119             ENDDO
120            ELSE
121             DO k=1,Nr
122              DO j=1-Oly,sNy+Oly
123               DO i=1-Olx,sNx+Olx
124                localTijk(i,j,k) = tracer(i,j,k,bi,bj)
125               ENDDO
126              ENDDO
127             ENDDO
128            ENDIF
129          ENDIF
130    
131        IF (implicitDiffusion) THEN        IF (implicitDiffusion) THEN
132  C--   set the tri-diagonal matrix to solve the implicit diffusion problem  C--   set the tri-diagonal matrix to solve the implicit diffusion problem
133         diagonalNumber = 3         diagonalNumber = 3
# Line 83  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 94  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 115  C--   end if implicitDiffusion Line 165  C--   end if implicitDiffusion
165    
166        IF (implicitAdvection) THEN        IF (implicitAdvection) THEN
167    
168         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  
169    
170            DO j=1-Oly,sNy+Oly  C--    Compute transport
171             DO i=1-Olx,sNx+Olx          IF (k.EQ.Nr) THEN
172              rTrans(i,j) = wVel(i,j,k,bi,bj)*rA(i,j,bi,bj)           DO j=1-Oly,sNy+Oly
173       &                   *maskC(i,j,k-1,bi,bj)            DO i=1-Olx,sNx+Olx
174             ENDDO              rTransKp1(i,j) = 0.
175              ENDDO
176             ENDDO
177            ELSE
178             DO j=1-Oly,sNy+Oly
179              DO i=1-Olx,sNx+Olx
180                rTransKp1(i,j) = rTrans(i,j)
181              ENDDO
182             ENDDO
183            ENDIF
184    
185            IF (k.EQ.1) THEN
186             DO j=1-Oly,sNy+Oly
187              DO i=1-Olx,sNx+Olx
188                wFld(i,j)   = 0. _d 0
189                rTrans(i,j) = 0. _d 0
190            ENDDO            ENDDO
191             ENDDO
192            ELSE
193             DO j=1-Oly,sNy+Oly
194              DO i=1-Olx,sNx+Olx
195                wFld(i,j)   = wVel(i,j,k,bi,bj)
196                rTrans(i,j) = wFld(i,j)*rA(i,j,bi,bj)*maskC(i,j,k-1,bi,bj)
197              ENDDO
198             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
207  C-       space Centered advection scheme, Flux form:          DO j=jMin,jMax
208            DO j=jMin,jMax            DO i=iMin,iMax
209             DO i=iMin,iMax  c          localTijk(i,j,k)     = gTracer(i,j,k,bi,bj)
210              rFlx = 0.5 _d 0 *deltaTtracer*rTrans(i,j)             gTracer(i,j,k,bi,bj) = gTracer(i,j,k,bi,bj)
211       &                      *recip_rA(i,j,bi,bj)*rkFac       &      + deltaTLev(1)*recip_rA(i,j,bi,bj)
212              b5d(i,j,k) = b5d(i,j,k)       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
213       &                 + rFlx*recip_hFacC(i,j,k,bi,bj)*recip_drF(k)       &       *tracer(i,j,k,bi,bj)*(rTransKp1(i,j)-rTrans(i,j))*rkSign
             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  
214            ENDDO            ENDDO
215            ENDDO
216    
217  C--      end k loop  #ifdef ALLOW_AIM
218           ENDDO  C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr
219         ELSE          IF ( K.GE.2 .AND.
220           STOP 'GAD_IMPLICIT_R: Adv.Scheme in Impl form not yet coded'       &     (.NOT.useAIM .OR.tracerIdentity.NE.GAD_SALINITY .OR.K.LT.Nr)
221         ENDIF       &              ) THEN
222    #else
223            IF ( K.GE.2 ) THEN
224    #endif
225    
226             IF ( advectionScheme.EQ.ENUM_CENTERED_2ND ) THEN
227              diagonalNumber = 3
228              CALL GAD_C2_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,
229         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,
238         I                        myThid )
239             ELSEIF ( advectionScheme.EQ.ENUM_FLUX_LIMIT ) THEN
240              diagonalNumber = 3
241              CALL GAD_FLUXLIMIT_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,
242         I                        deltaTLev, rTrans, localTijk,
243         U                        b5d, c5d, d5d,
244         I                        myThid )
245             ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_3RD
246         &       .OR. advectionScheme.EQ.ENUM_CENTERED_4TH
247         &       .OR. advectionScheme.EQ.ENUM_DST3 ) THEN
248              diagonalNumber = 5
249              CALL GAD_U3C4_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,
250         I                        advectionScheme, deltaTLev, rTrans,
251         U                        a5d, b5d, c5d, d5d, e5d,
252         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
260              STOP 'GAD_IMPLICIT_R: Adv.Scheme in Impl form not yet coded'
261             ENDIF
262    
263            ENDIF
264    
265    C--     end k loop
266           ENDDO
267    
268  C--   end if implicitAdvection  C--   end if implicitAdvection
269        ENDIF        ENDIF
# Line 168  C--   Solve tri-diagonal system : Line 278  C--   Solve tri-diagonal system :
278          IF (errCode.GE.1) THEN          IF (errCode.GE.1) THEN
279            STOP 'GAD_IMPLICIT_R: error when solving 3-Diag problem'            STOP 'GAD_IMPLICIT_R: error when solving 3-Diag problem'
280          ENDIF          ENDIF
281          ELSEIF ( diagonalNumber .EQ. 5 ) THEN
282    C--   Solve penta-diagonal system :
283            CALL SOLVE_PENTADIAGONAL( iMin,iMax, jMin,jMax,
284         I                            a5d, b5d, c5d, d5d, e5d,
285         U                            gTracer,
286         O                            errCode,
287         I                            bi, bj, myThid )
288            IF (errCode.GE.1) THEN
289              STOP 'GAD_IMPLICIT_R: error when solving 5-Diag problem'
290            ENDIF
291        ELSEIF ( diagonalNumber .NE. 1 ) THEN        ELSEIF ( diagonalNumber .NE. 1 ) THEN
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.1  
changed lines
  Added in v.1.14
  ViewVC Help
Powered by ViewVC 1.1.22