pkg/generic_advdiff/gad_implicit_r.F

C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_implicit_r.F,v 1.1 2004/01/03 00:43:01 jmc Exp $
C $Name:  $

#include "GAD_OPTIONS.h"

CBOP
C     !ROUTINE: GAD_IMPLICIT_R
C     !INTERFACE:
      SUBROUTINE GAD_IMPLICIT_R( 
     I      implicitAdvection, advectionScheme, tracerIdentity,
     I      kappaRX, wVel, tracer, 
     U      gTracer,
     I      bi, bj, myTime, myIter, myThid )
C     !DESCRIPTION: \bv
C     *==========================================================*
C     | S/R GAD_IMPLICIT_R
C     | o Solve implicitly vertical advection & diffusion
C     |   for one tracer
C     *==========================================================*
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE
C     == Global data ==
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "GAD.h"

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine Arguments ==
      LOGICAL implicitAdvection
      INTEGER advectionScheme
      INTEGER tracerIdentity
      _RL kappaRX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
      _RL wVel   (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
      _RL tracer (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
      _RL gTracer(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
      INTEGER bi, bj
      _RL myTime
      INTEGER myIter, myThid

C     !LOCAL VARIABLES:
C     == Local variables ==
      INTEGER iMin,iMax,jMin,jMax
      INTEGER i,j,k
      INTEGER diagonalNumber, errCode
      _RL a5d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
      _RL b5d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
      _RL c5d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
      _RL d5d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
      _RL e5d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
      _RL rTrans(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      _RL rTransKp1(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      _RL localTijk(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
      _RL limitDf  (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      _RL rFlx
CEOP

      IF (Nr.LE.1) RETURN

C--   Initialise
      iMin = 1
      jMin = 1
      iMax = sNx
      jMax = sNy
      DO k=1,Nr
       DO j=1-Oly,sNy+Oly
        DO i=1-Olx,sNx+Olx
         a5d(i,j,k) = 0. _d 0
         b5d(i,j,k) = 0. _d 0
         c5d(i,j,k) = 1. _d 0
         d5d(i,j,k) = 0. _d 0
         e5d(i,j,k) = 0. _d 0
        ENDDO
       ENDDO
      ENDDO
      diagonalNumber = 1

C--   Non-Linear Advection scheme: keep a local copy of tracer field
      IF ( advectionScheme.EQ.ENUM_FLUX_LIMIT .OR. 
     &     advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
        IF ( multiDimAdvection ) THEN
         DO k=1,Nr
          DO j=1-Oly,sNy+Oly
           DO i=1-Olx,sNx+Olx
            localTijk(i,j,k) = gTracer(i,j,k,bi,bj)
           ENDDO
          ENDDO
         ENDDO
        ELSE
         DO k=1,Nr
          DO j=1-Oly,sNy+Oly
           DO i=1-Olx,sNx+Olx
            localTijk(i,j,k) = tracer(i,j,k,bi,bj)
           ENDDO
          ENDDO
         ENDDO
        ENDIF
      ENDIF

      IF (implicitDiffusion) THEN
C--   set the tri-diagonal matrix to solve the implicit diffusion problem
       diagonalNumber = 3
C-     1rst lower diagonal :
       DO k=2,Nr
        DO j=jMin,jMax
         DO i=iMin,iMax
          IF (maskC(i,j,k-1,bi,bj).EQ.1.) 
     &     b5d(i,j,k) = -deltaTtracer
     &                  *recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
     &                  *kappaRX(i,j, k )*recip_drC( k )
         ENDDO
        ENDDO
       ENDDO
C-     1rst upper diagonal :
       DO k=1,Nr-1
        DO j=jMin,jMax
         DO i=iMin,iMax
          IF (maskC(i,j,k+1,bi,bj).EQ.1.) 
     &     d5d(i,j,k) = -deltaTtracer
     &                 *recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
     &                 *KappaRX(i,j,k+1)*recip_drC(k+1)
         ENDDO
        ENDDO
       ENDDO
C-     Main diagonal :
       DO k=1,Nr
        DO j=jMin,jMax
         DO i=iMin,iMax
           c5d(i,j,k) = 1. _d 0 - b5d(i,j,k) - d5d(i,j,k)
         ENDDO
        ENDDO
       ENDDO

C--   end if implicitDiffusion
      ENDIF

      IF (implicitAdvection) THEN

       DO k=Nr,1,-1

C--    Compute transport
        IF (k.EQ.Nr) THEN
         DO j=1-Oly,sNy+Oly
          DO i=1-Olx,sNx+Olx
            rTransKp1(i,j) = 0.
          ENDDO
         ENDDO
        ELSE
         DO j=1-Oly,sNy+Oly
          DO i=1-Olx,sNx+Olx
            rTransKp1(i,j) = rTrans(i,j)
          ENDDO
         ENDDO
        ENDIF

        IF (k.EQ.1) THEN
         DO j=1-Oly,sNy+Oly
          DO i=1-Olx,sNx+Olx
            rTrans(i,j) = 0.
          ENDDO
         ENDDO
        ELSE
         DO j=1-Oly,sNy+Oly
          DO i=1-Olx,sNx+Olx
            rTrans(i,j) = wVel(i,j,k,bi,bj)*rA(i,j,bi,bj)
     &                   *maskC(i,j,k-1,bi,bj)
          ENDDO
         ENDDO
#ifdef ALLOW_GMREDI
C--   Residual transp = Bolus transp + Eulerian transp
          IF (useGMRedi)
     &      CALL GMREDI_CALC_WFLOW(
     &                    rTrans, bi, bj, k, myThid)
#endif /* ALLOW_GMREDI */
        ENDIF
        DO j=jMin,jMax
          DO i=iMin,iMax
c          localTijk(i,j,k)     = gTracer(i,j,k,bi,bj) 
           gTracer(i,j,k,bi,bj) = gTracer(i,j,k,bi,bj)
     &      + deltaTtracer*recip_rA(i,j,bi,bj)
     &       *_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
     &       *tracer(i,j,k,bi,bj)*(rTrans(i,j)-rTransKp1(i,j))*rkFac
          ENDDO
        ENDDO

        IF (K.GE.2) THEN

         IF ( advectionScheme.EQ.ENUM_CENTERED_2ND ) THEN
          diagonalNumber = 3
          CALL GAD_C2_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,
     I                        deltaTtracer, rTrans,
     U                        b5d, c5d, d5d,
     I                        myThid)
         ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
          diagonalNumber = 3
          CALL GAD_FLUXLIMIT_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,
     I                        deltaTtracer, rTrans, localTijk,
     U                        b5d, c5d, d5d,
     I                        myThid)
         ELSEIF (advectionScheme.EQ.ENUM_UPWIND_3RD .OR.
     &           advectionScheme.EQ.ENUM_CENTERED_4TH) THEN
          diagonalNumber = 5
          CALL GAD_U3C4_IMPL_R( bi,bj,k, iMin,iMax,jMin,jMax,
     I                        advectionScheme, deltaTtracer, rTrans,
     U                        a5d, b5d, c5d, d5d, e5d,
     I                        myThid)
         ELSE
          STOP 'GAD_IMPLICIT_R: Adv.Scheme in Impl form not yet coded'
         ENDIF

        ENDIF

C--     end k loop
       ENDDO

C--   end if implicitAdvection
      ENDIF

      IF ( diagonalNumber .EQ. 3 ) THEN
C--   Solve tri-diagonal system :
        CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,
     I                          b5d, c5d, d5d,
     U                          gTracer,
     O                          errCode,
     I                          bi, bj, myThid )
        IF (errCode.GE.1) THEN
          STOP 'GAD_IMPLICIT_R: error when solving 3-Diag problem'
        ENDIF
      ELSEIF ( diagonalNumber .EQ. 5 ) THEN
C--   Solve penta-diagonal system :
        CALL SOLVE_PENTADIAGONAL( iMin,iMax, jMin,jMax,
     I                            a5d, b5d, c5d, d5d, e5d,
     U                            gTracer,
     O                            errCode,
     I                            bi, bj, myThid )
        IF (errCode.GE.1) THEN
          STOP 'GAD_IMPLICIT_R: error when solving 5-Diag problem'
        ENDIF
      ELSEIF ( diagonalNumber .NE. 1 ) THEN
        STOP 'GAD_IMPLICIT_R: no solver available'
      ENDIF

      RETURN
      END
1	jmc	1.2	C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_implicit_r.F,v 1.1 2004/01/03 00:43:01 jmc Exp $
2	jmc	1.1	C $Name: $
3
4			#include "GAD_OPTIONS.h"
5
6			CBOP
7			C !ROUTINE: GAD_IMPLICIT_R
8			C !INTERFACE:
9			SUBROUTINE GAD_IMPLICIT_R(
10			I implicitAdvection, advectionScheme, tracerIdentity,
11			I kappaRX, wVel, tracer,
12			U gTracer,
13			I bi, bj, myTime, myIter, myThid )
14			C !DESCRIPTION: \bv
15			C ==========================================================
16			C \| S/R GAD_IMPLICIT_R
17			C \| o Solve implicitly vertical advection & diffusion
18			C \| for one tracer
19			C ==========================================================
20			C ==========================================================
21			C \ev
22
23			C !USES:
24			IMPLICIT NONE
25			C == Global data ==
26			#include "SIZE.h"
27			#include "EEPARAMS.h"
28			#include "PARAMS.h"
29			#include "GRID.h"
30			#include "GAD.h"
31
32			C !INPUT/OUTPUT PARAMETERS:
33			C == Routine Arguments ==
34			LOGICAL implicitAdvection
35			INTEGER advectionScheme
36			INTEGER tracerIdentity
37			_RL kappaRX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
38			_RL wVel (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
39			_RL tracer (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
40			_RL gTracer(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
41			INTEGER bi, bj
42			_RL myTime
43			INTEGER myIter, myThid
44
45			C !LOCAL VARIABLES:
46			C == Local variables ==
47			INTEGER iMin,iMax,jMin,jMax
48			INTEGER i,j,k
49			INTEGER diagonalNumber, errCode
50			_RL a5d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
51			_RL b5d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
52			_RL c5d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
53			_RL d5d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
54			_RL e5d(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
55			_RL rTrans(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
56	jmc	1.2	_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	jmc	1.1	_RL rFlx
60			CEOP
61
62			IF (Nr.LE.1) RETURN
63
64			C-- Initialise
65			iMin = 1
66			jMin = 1
67			iMax = sNx
68			jMax = sNy
69			DO k=1,Nr
70			DO j=1-Oly,sNy+Oly
71			DO i=1-Olx,sNx+Olx
72			a5d(i,j,k) = 0. _d 0
73			b5d(i,j,k) = 0. _d 0
74			c5d(i,j,k) = 1. _d 0
75			d5d(i,j,k) = 0. _d 0
76			e5d(i,j,k) = 0. _d 0
77			ENDDO
78			ENDDO
79			ENDDO
80			diagonalNumber = 1
81
82	jmc	1.2	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	jmc	1.1	IF (implicitDiffusion) THEN
105			C-- set the tri-diagonal matrix to solve the implicit diffusion problem
106			diagonalNumber = 3
107			C- 1rst lower diagonal :
108			DO k=2,Nr
109			DO j=jMin,jMax
110			DO i=iMin,iMax
111			IF (maskC(i,j,k-1,bi,bj).EQ.1.)
112			& b5d(i,j,k) = -deltaTtracer
113			& recip_hFacC(i,j,k,bi,bj)recip_drF(k)
114			& kappaRX(i,j, k )recip_drC( k )
115			ENDDO
116			ENDDO
117			ENDDO
118			C- 1rst upper diagonal :
119			DO k=1,Nr-1
120			DO j=jMin,jMax
121			DO i=iMin,iMax
122			IF (maskC(i,j,k+1,bi,bj).EQ.1.)
123			& d5d(i,j,k) = -deltaTtracer
124			& recip_hFacC(i,j,k,bi,bj)recip_drF(k)
125			& KappaRX(i,j,k+1)recip_drC(k+1)
126			ENDDO
127			ENDDO
128			ENDDO
129			C- Main diagonal :
130			DO k=1,Nr
131			DO j=jMin,jMax
132			DO i=iMin,iMax
133			c5d(i,j,k) = 1. _d 0 - b5d(i,j,k) - d5d(i,j,k)
134			ENDDO
135			ENDDO
136			ENDDO
137
138			C-- end if implicitDiffusion
139			ENDIF
140
141			IF (implicitAdvection) THEN
142
143	jmc	1.2	DO k=Nr,1,-1
144
145			C-- Compute transport
146			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	jmc	1.1
160	jmc	1.2	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	jmc	1.1	rTrans(i,j) = wVel(i,j,k,bi,bj)*rA(i,j,bi,bj)
170			& *maskC(i,j,k-1,bi,bj)
171			ENDDO
172	jmc	1.2	ENDDO
173	jmc	1.1	#ifdef ALLOW_GMREDI
174			C-- Residual transp = Bolus transp + Eulerian transp
175			IF (useGMRedi)
176			& CALL GMREDI_CALC_WFLOW(
177			& rTrans, bi, bj, k, myThid)
178			#endif /* ALLOW_GMREDI */
179	jmc	1.2	ENDIF
180			DO j=jMin,jMax
181			DO i=iMin,iMax
182			c localTijk(i,j,k) = gTracer(i,j,k,bi,bj)
183			gTracer(i,j,k,bi,bj) = gTracer(i,j,k,bi,bj)
184			& + deltaTtracer*recip_rA(i,j,bi,bj)
185			& _recip_hFacC(i,j,k,bi,bj)recip_drF(k)
186			& tracer(i,j,k,bi,bj)(rTrans(i,j)-rTransKp1(i,j))*rkFac
187			ENDDO
188			ENDDO
189
190			IF (K.GE.2) THEN
191	jmc	1.1
192	jmc	1.2	IF ( advectionScheme.EQ.ENUM_CENTERED_2ND ) THEN
193			diagonalNumber = 3
194			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	jmc	1.1
217	jmc	1.2	C-- end k loop
218			ENDDO
219	jmc	1.1
220			C-- end if implicitAdvection
221			ENDIF
222
223			IF ( diagonalNumber .EQ. 3 ) THEN
224			C-- Solve tri-diagonal system :
225			CALL SOLVE_TRIDIAGONAL( iMin,iMax, jMin,jMax,
226			I b5d, c5d, d5d,
227			U gTracer,
228			O errCode,
229			I bi, bj, myThid )
230			IF (errCode.GE.1) THEN
231			STOP 'GAD_IMPLICIT_R: error when solving 3-Diag problem'
232	jmc	1.2	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	jmc	1.1	ENDIF
243			ELSEIF ( diagonalNumber .NE. 1 ) THEN
244			STOP 'GAD_IMPLICIT_R: no solver available'
245			ENDIF
246
247			RETURN
248			END