/[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.12 by heimbach, Wed Jun 7 01:55:14 2006 UTC
# Line 11  C     !INTERFACE: Line 11  C     !INTERFACE:
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: \bv  C     !DESCRIPTION:
15  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  
16    
17  C     !USES:  C     !USES:
18        IMPLICIT NONE        IMPLICIT NONE
# Line 29  C     == Global data == Line 23  C     == Global data ==
23  #include "GRID.h"  #include "GRID.h"
24  #include "GAD.h"  #include "GAD.h"
25    
26  C     !INPUT/OUTPUT PARAMETERS:  C !INPUT/OUTPUT PARAMETERS:
27  C     == Routine Arguments ==  C == Routine Arguments ==
28    C implicitAdvection :: if True, treat vertical advection implicitly
29    C advectionScheme   :: advection scheme to use
30    C tracerIdentity    :: Identifier for the tracer
31    C kappaRX           :: 3-D array for vertical diffusion coefficient
32    C wVel              :: vertical component of the velcity field
33    C tracer            :: tracer field at current time step
34    C gTracer           :: future tracer field
35    C bi,bj             :: tile indices
36    C myTime            :: current time
37    C myIter            :: current iteration number
38    C myThid            :: thread number
39        LOGICAL implicitAdvection        LOGICAL implicitAdvection
40        INTEGER advectionScheme        INTEGER advectionScheme
41        INTEGER tracerIdentity        INTEGER tracerIdentity
# Line 42  C     == Routine Arguments == Line 47  C     == Routine Arguments ==
47        _RL myTime        _RL myTime
48        INTEGER myIter, myThid        INTEGER myIter, myThid
49    
50  C     !LOCAL VARIABLES:  C !LOCAL VARIABLES:
51  C     == Local variables ==  C == Local variables ==
52    C iMin,iMax,jMin,jMax  :: computational domain
53    C i,j,k  :: loop indices
54    C a5d    :: 2nd  lower diagonal of the pentadiagonal matrix
55    C b5d    :: 1rst lower diagonal of the pentadiagonal matrix
56    C c5d    :: main diagonal       of the pentadiagonal matrix
57    C d5d    :: 1rst upper diagonal of the pentadiagonal matrix
58    C e5d    :: 2nd  upper diagonal of the pentadiagonal matrix
59    C rTrans    :: vertical volume transport at inteface k
60    C rTransKp1 :: vertical volume transport at inteface k+1
61    C localTijk :: local copy of tracer (for Non-Lin Adv.Scheme)
62    C diagonalNumber :: number of non-zero diagonals in the matrix
63    C errCode   :: > 0 if singular matrix
64        INTEGER iMin,iMax,jMin,jMax        INTEGER iMin,iMax,jMin,jMax
65        INTEGER i,j,k        INTEGER i,j,k
66        INTEGER diagonalNumber, errCode        INTEGER diagonalNumber, errCode
# Line 53  C     == Local variables == Line 70  C     == Local variables ==
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 rTrans(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
73        _RL rFlx        _RL rTransKp1(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
74          _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 76  C--   Initialise Line 103  C--   Initialise
103        ENDDO        ENDDO
104        diagonalNumber = 1        diagonalNumber = 1
105    
106    C--   Non-Linear Advection scheme: keep a local copy of tracer field
107          IF ( advectionScheme.EQ.ENUM_FLUX_LIMIT .OR.
108         &     advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
109            IF ( multiDimAdvection ) THEN
110             DO k=1,Nr
111              DO j=1-Oly,sNy+Oly
112               DO i=1-Olx,sNx+Olx
113                localTijk(i,j,k) = gTracer(i,j,k,bi,bj)
114               ENDDO
115              ENDDO
116             ENDDO
117            ELSE
118             DO k=1,Nr
119              DO j=1-Oly,sNy+Oly
120               DO i=1-Olx,sNx+Olx
121                localTijk(i,j,k) = tracer(i,j,k,bi,bj)
122               ENDDO
123              ENDDO
124             ENDDO
125            ENDIF
126          ENDIF
127    
128        IF (implicitDiffusion) THEN        IF (implicitDiffusion) THEN
129  C--   set the tri-diagonal matrix to solve the implicit diffusion problem  C--   set the tri-diagonal matrix to solve the implicit diffusion problem
130         diagonalNumber = 3         diagonalNumber = 3
# Line 83  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)
136       &     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)  
137       &                  *kappaRX(i,j, k )*recip_drC( k )       &                  *kappaRX(i,j, k )*recip_drC( k )
138           ENDDO           ENDDO
139          ENDDO          ENDDO
# Line 94  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)
146       &     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)  
147       &                 *KappaRX(i,j,k+1)*recip_drC(k+1)       &                 *KappaRX(i,j,k+1)*recip_drC(k+1)
148           ENDDO           ENDDO
149          ENDDO          ENDDO
# Line 115  C--   end if implicitDiffusion Line 162  C--   end if implicitDiffusion
162    
163        IF (implicitAdvection) THEN        IF (implicitAdvection) THEN
164    
165         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  
166    
167            DO j=1-Oly,sNy+Oly  C--    Compute transport
168             DO i=1-Olx,sNx+Olx          IF (k.EQ.Nr) THEN
169             DO j=1-Oly,sNy+Oly
170              DO i=1-Olx,sNx+Olx
171                rTransKp1(i,j) = 0.
172              ENDDO
173             ENDDO
174            ELSE
175             DO j=1-Oly,sNy+Oly
176              DO i=1-Olx,sNx+Olx
177                rTransKp1(i,j) = rTrans(i,j)
178              ENDDO
179             ENDDO
180            ENDIF
181    
182            IF (k.EQ.1) THEN
183             DO j=1-Oly,sNy+Oly
184              DO i=1-Olx,sNx+Olx
185                rTrans(i,j) = 0.
186              ENDDO
187             ENDDO
188            ELSE
189             DO j=1-Oly,sNy+Oly
190              DO i=1-Olx,sNx+Olx
191              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)
192       &                   *maskC(i,j,k-1,bi,bj)       &                   *maskC(i,j,k-1,bi,bj)
            ENDDO  
193            ENDDO            ENDDO
194             ENDDO
195  #ifdef ALLOW_GMREDI  #ifdef ALLOW_GMREDI
196  C--   Residual transp = Bolus transp + Eulerian transp  C--   Residual transp = Bolus transp + Eulerian transp
197            IF (useGMRedi)            IF (useGMRedi)
198       &      CALL GMREDI_CALC_WFLOW(       &      CALL GMREDI_CALC_WFLOW(
199       &                    rTrans, bi, bj, k, myThid)       &                    rTrans, bi, bj, k, myThid)
200  #endif /* ALLOW_GMREDI */  #endif /* ALLOW_GMREDI */
201            ENDIF
202  C-       space Centered advection scheme, Flux form:          DO j=jMin,jMax
203            DO j=jMin,jMax            DO i=iMin,iMax
204             DO i=iMin,iMax  c          localTijk(i,j,k)     = gTracer(i,j,k,bi,bj)
205              rFlx = 0.5 _d 0 *deltaTtracer*rTrans(i,j)             gTracer(i,j,k,bi,bj) = gTracer(i,j,k,bi,bj)
206       &                      *recip_rA(i,j,bi,bj)*rkFac       &      + dTtracerLev(1)*recip_rA(i,j,bi,bj)
207              b5d(i,j,k) = b5d(i,j,k)       &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
208       &                 + 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  
209            ENDDO            ENDDO
210            ENDDO
211    
212  C--      end k loop  #ifdef ALLOW_AIM
213           ENDDO  C- a hack to prevent Water-Vapor vert.transport into the stratospheric level Nr
214         ELSE          IF ( K.GE.2 .AND.
215           STOP 'GAD_IMPLICIT_R: Adv.Scheme in Impl form not yet coded'       &     (.NOT.useAIM .OR.tracerIdentity.NE.GAD_SALINITY .OR.K.LT.Nr)
216         ENDIF       &              ) THEN
217    #else
218            IF ( K.GE.2 ) THEN
219    #endif
220    
221             IF ( advectionScheme.EQ.ENUM_CENTERED_2ND ) THEN
222              diagonalNumber = 3
223              CALL GAD_C2_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,
224         I                        dTtracerLev, rTrans,
225         U                        b5d, c5d, d5d,
226         I                        myThid )
227             ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_1RST
228         &       .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
229              diagonalNumber = 3
230              CALL GAD_DST2U1_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,
231         I                        advectionScheme, dTtracerLev, rTrans,
232         U                        b5d, c5d, d5d,
233         I                        myThid )
234             ELSEIF ( advectionScheme.EQ.ENUM_FLUX_LIMIT ) THEN
235              diagonalNumber = 3
236              CALL GAD_FLUXLIMIT_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,
237         I                        dTtracerLev, rTrans, localTijk,
238         U                        b5d, c5d, d5d,
239         I                        myThid )
240             ELSEIF ( advectionScheme.EQ.ENUM_UPWIND_3RD
241         &       .OR. advectionScheme.EQ.ENUM_CENTERED_4TH
242         &       .OR. advectionScheme.EQ.ENUM_DST3 ) THEN
243              diagonalNumber = 5
244              CALL GAD_U3C4_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,
245         I                        advectionScheme, dTtracerLev, rTrans,
246         U                        a5d, b5d, c5d, d5d, e5d,
247         I                        myThid )
248             ELSEIF ( advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
249              diagonalNumber = 5
250              CALL GAD_DST3FL_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,
251         I                        dTtracerLev, rTrans, localTijk,
252         U                        a5d, b5d, c5d, d5d, e5d,
253         I                        myThid )
254             ELSE
255              STOP 'GAD_IMPLICIT_R: Adv.Scheme in Impl form not yet coded'
256             ENDIF
257    
258            ENDIF
259    
260    C--     end k loop
261           ENDDO
262    
263  C--   end if implicitAdvection  C--   end if implicitAdvection
264        ENDIF        ENDIF
# Line 168  C--   Solve tri-diagonal system : Line 273  C--   Solve tri-diagonal system :
273          IF (errCode.GE.1) THEN          IF (errCode.GE.1) THEN
274            STOP 'GAD_IMPLICIT_R: error when solving 3-Diag problem'            STOP 'GAD_IMPLICIT_R: error when solving 3-Diag problem'
275          ENDIF          ENDIF
276          ELSEIF ( diagonalNumber .EQ. 5 ) THEN
277    C--   Solve penta-diagonal system :
278            CALL SOLVE_PENTADIAGONAL( iMin,iMax, jMin,jMax,
279         I                            a5d, b5d, c5d, d5d, e5d,
280         U                            gTracer,
281         O                            errCode,
282         I                            bi, bj, myThid )
283            IF (errCode.GE.1) THEN
284              STOP 'GAD_IMPLICIT_R: error when solving 5-Diag problem'
285            ENDIF
286        ELSEIF ( diagonalNumber .NE. 1 ) THEN        ELSEIF ( diagonalNumber .NE. 1 ) THEN
287          STOP 'GAD_IMPLICIT_R: no solver available'          STOP 'GAD_IMPLICIT_R: no solver available'
288        ENDIF        ENDIF
289    
290    #ifdef ALLOW_DIAGNOSTICS
291    C--   Set diagnostic suffix for the current tracer
292          IF ( useDiagnostics .AND. implicitDiffusion ) THEN
293            diagSufx = GAD_DIAG_SUFX( tracerIdentity, myThid )
294            diagName = 'DFrI'//diagSufx
295            IF ( DIAGNOSTICS_IS_ON(diagName,myThid) ) THEN
296             DO k= 1,Nr
297              IF ( k.EQ.1 ) THEN
298    C-  Note: Needs to call DIAGNOSTICS_FILL at level k=1 even if array == 0
299    C         otherwise counter is not incremented !!
300                DO j=1-OLy,sNy+OLy
301                 DO i=1-OLx,sNx+OLx
302                   df(i,j) = 0. _d 0
303                 ENDDO
304                ENDDO
305              ELSE
306                DO j=1,sNy
307                 DO i=1,sNx
308                   df(i,j) =
309         &              rA(i,j,bi,bj)
310         &            * KappaRX(i,j,k)*recip_drC(k)
311         &            * (gTracer(i,j,k,bi,bj) - gTracer(i,j,k-1,bi,bj))
312                 ENDDO
313                ENDDO
314              ENDIF
315              CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
316             ENDDO
317            ENDIF
318          ENDIF
319    #endif /* ALLOW_DIAGNOSTICS */
320    
321    C--   end if Nr > 1
322          ENDIF
323    
324        RETURN        RETURN
325        END        END
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.12
  ViewVC Help
Powered by ViewVC 1.1.22