model/src/impldiff.F

C $Header: /u/gcmpack/MITgcm/model/src/impldiff.F,v 1.27 2009/06/26 23:10:09 jahn Exp $
C $Name:  $

#include "PACKAGES_CONFIG.h"
#include "CPP_OPTIONS.h"

CBOP
C     !ROUTINE: IMPLDIFF
C     !INTERFACE:
      SUBROUTINE IMPLDIFF( bi, bj, iMin, iMax, jMin, jMax,
     I                     tracerId, KappaRX, recip_hFac,
     U                     gXNm1,
     I                     myThid )
C     !DESCRIPTION: \bv
C     *==========================================================*
C     | S/R IMPLDIFF                                              
C     | o Solve implicit diffusion equation for vertical          
C     |   diffusivity.                                            
C     *==========================================================*
C     | o Recoded from 2d intermediate fields to 3d to reduce     
C     |   TAMC storage                                            
C     | o Fixed missing masks for fields a(), c()                 
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE
C     == Global data ==
#include "SIZE.h"
#include "DYNVARS.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#ifdef ALLOW_GENERIC_ADVDIFF
#include "GAD.h"
#endif
#ifdef ALLOW_LONGSTEP
#include "LONGSTEP_PARAMS.h"
#endif
#ifdef ALLOW_PTRACERS
#include "PTRACERS_SIZE.h"
#include "PTRACERS_PARAMS.h"
#endif
#ifdef ALLOW_AUTODIFF_TAMC
#include "tamc_keys.h"
#endif

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine Arguments ==
C     tracerId :: tracer Identificator (if > 0) ; 
C      = 0 or < 0 when solving vertical viscosity implicitly for U or V
      INTEGER bi,bj,iMin,iMax,jMin,jMax
      INTEGER tracerId
      _RL KappaRX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
      _RS recip_hFac(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
      _RL gXnm1(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
      INTEGER myThid

C     !LOCAL VARIABLES:
C     == Local variables ==
      INTEGER i,j,k
      _RL deltaTX(Nr)
      _RL gYnm1(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
      _RL a(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
      _RL b(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
      _RL c(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
      _RL bet(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
      _RL gam(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
#ifdef ALLOW_DIAGNOSTICS
      CHARACTER*8 diagName
      CHARACTER*4 diagSufx
#ifdef ALLOW_GENERIC_ADVDIFF
      CHARACTER*4 GAD_DIAG_SUFX
      EXTERNAL    GAD_DIAG_SUFX
#endif
      LOGICAL     DIAGNOSTICS_IS_ON
      EXTERNAL    DIAGNOSTICS_IS_ON
      _RL df (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
#endif /* ALLOW_DIAGNOSTICS */
CEOP

cph(
cph Not good for TAF: may create irreducible control flow graph
cph      IF (Nr.LE.1) RETURN
cph)

#ifdef ALLOW_PTRACERS
      IF ( tracerId.GE.GAD_TR1) THEN
        DO k=1,Nr
         deltaTX(k) = PTRACERS_dTLev(k)
        ENDDO
      ELSEIF ( tracerId.GE.1 ) THEN
#else
      IF ( tracerId.GE.1 ) THEN
#endif
        DO k=1,Nr
         deltaTX(k) = dTtracerLev(k)
        ENDDO
      ELSE
        DO k=1,Nr
         deltaTX(k) = deltaTmom
        ENDDO
      ENDIF

C--   Initialise
      DO k=1,Nr
#ifdef TARGET_NEC_SX
       DO j=1-OLy,sNy+OLy
        DO i=1-OLx,sNx+OLx
#else
       DO j=jMin,jMax
        DO i=iMin,iMax
#endif
         gYNm1(i,j,k,bi,bj) = 0. _d 0
        ENDDO
       ENDDO
      ENDDO

C--   Old aLower
#ifdef TARGET_NEC_SX
      DO j=1-OLy,sNy+OLy
       DO i=1-OLx,sNx+OLx
#else
      DO j=jMin,jMax
       DO i=iMin,iMax
#endif
         a(i,j,1) = 0. _d 0 
       ENDDO
      ENDDO
      DO k=2,Nr
#ifdef TARGET_NEC_SX
       DO j=1-OLy,sNy+OLy
        DO i=1-OLx,sNx+OLx
#else
       DO j=jMin,jMax
        DO i=iMin,iMax
#endif
          a(i,j,k) = -deltaTX(k)*recip_hFac(i,j,k,bi,bj)*recip_drF(k)
     &               *recip_deepFac2C(k)*recip_rhoFacC(k)
     &               *KappaRX(i,j, k )*recip_drC( k )
     &               *deepFac2F(k)*rhoFacF(k)
          IF (recip_hFac(i,j,k-1,bi,bj).EQ.0.) a(i,j,k)=0.
        ENDDO
       ENDDO
      ENDDO

C--   Old aUpper
      DO k=1,Nr-1
#ifdef TARGET_NEC_SX
       DO j=1-OLy,sNy+OLy
        DO i=1-OLx,sNx+OLx
#else
       DO j=jMin,jMax
        DO i=iMin,iMax
#endif
          c(i,j,k) = -deltaTX(k)*recip_hFac(i,j,k,bi,bj)*recip_drF(k)
     &               *recip_deepFac2C(k)*recip_rhoFacC(k)
     &               *KappaRX(i,j,k+1)*recip_drC(k+1)
     &               *deepFac2F(k+1)*rhoFacF(k+1)
          IF (recip_hFac(i,j,k+1,bi,bj).EQ.0.) c(i,j,k)=0.
        ENDDO
       ENDDO
      ENDDO
#ifdef TARGET_NEC_SX
      DO j=1-OLy,sNy+OLy
       DO i=1-OLx,sNx+OLx
#else
      DO j=jMin,jMax
       DO i=iMin,iMax
#endif
         c(i,j,Nr) = 0. _d 0
       ENDDO
      ENDDO

C--   Old aCenter
      DO k=1,Nr
#ifdef TARGET_NEC_SX
       DO j=1-OLy,sNy+OLy
        DO i=1-OLx,sNx+OLx
#else
       DO j=jMin,jMax
        DO i=iMin,iMax
#endif
          b(i,j,k) = 1. _d 0 - c(i,j,k) - a(i,j,k)
        ENDDO
       ENDDO
      ENDDO

C--   Old and new gam, bet are the same
      DO k=1,Nr
#ifdef TARGET_NEC_SX
       DO j=1-OLy,sNy+OLy
        DO i=1-OLx,sNx+OLx
#else
       DO j=jMin,jMax
        DO i=iMin,iMax
#endif
          bet(i,j,k) = 1. _d 0
          gam(i,j,k) = 0. _d 0
        ENDDO
       ENDDO
      ENDDO

C--   Only need do anything if Nr>1
      IF (Nr.GT.1) THEN

       k = 1
C--    Beginning of forward sweep (top level)
#ifdef TARGET_NEC_SX
       DO j=1-OLy,sNy+OLy
        DO i=1-OLx,sNx+OLx
#else
       DO j=jMin,jMax
        DO i=iMin,iMax
#endif
         IF (b(i,j,1).NE.0.) bet(i,j,1) = 1. _d 0 / b(i,j,1)
        ENDDO
       ENDDO

      ENDIF

C--   Middle of forward sweep
      IF (Nr.GE.2) THEN

CADJ loop = sequential
       DO k=2,Nr

#ifdef TARGET_NEC_SX
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
#else
        DO j=jMin,jMax
         DO i=iMin,iMax
#endif
          gam(i,j,k) = c(i,j,k-1)*bet(i,j,k-1)
          IF ( ( b(i,j,k) - a(i,j,k)*gam(i,j,k) ) .NE. 0.) 
     &        bet(i,j,k) = 1. _d 0 / ( b(i,j,k) - a(i,j,k)*gam(i,j,k) )
         ENDDO
        ENDDO

       ENDDO

      ENDIF


#ifdef TARGET_NEC_SX
      DO j=1-OLy,sNy+OLy
       DO i=1-OLx,sNx+OLx
#else
      DO j=jMin,jMax
       DO i=iMin,iMax
#endif
        gYNm1(i,j,1,bi,bj) = gXNm1(i,j,1,bi,bj)*bet(i,j,1)
       ENDDO
      ENDDO
      DO k=2,Nr
#ifdef TARGET_NEC_SX
       DO j=1-OLy,sNy+OLy
        DO i=1-OLx,sNx+OLx
#else
       DO j=jMin,jMax
        DO i=iMin,iMax
#endif
         gYnm1(i,j,k,bi,bj) = bet(i,j,k)*
     &        (gXnm1(i,j,k,bi,bj) - a(i,j,k)*gYnm1(i,j,k-1,bi,bj))
        ENDDO
       ENDDO
      ENDDO


C--    Backward sweep
CADJ loop = sequential
       DO k=Nr-1,1,-1
#ifdef TARGET_NEC_SX
            DO j=1-OLy,sNy+OLy
             DO i=1-OLx,sNx+OLx
#else
        DO j=jMin,jMax
         DO i=iMin,iMax
#endif
          gYnm1(i,j,k,bi,bj)=gYnm1(i,j,k,bi,bj)
     &              -gam(i,j,k+1)*gYnm1(i,j,k+1,bi,bj)
         ENDDO
        ENDDO
       ENDDO

       DO k=1,Nr
#ifdef TARGET_NEC_SX
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
#else
        DO j=jMin,jMax
         DO i=iMin,iMax
#endif
          gXnm1(i,j,k,bi,bj)=gYnm1(i,j,k,bi,bj)
         ENDDO
        ENDDO
       ENDDO

#ifdef ALLOW_DIAGNOSTICS
      IF ( useDiagnostics .AND.tracerId.NE.0 ) THEN
        IF ( tracerId.GE. 1 ) THEN
C--   Set diagnostic suffix for the current tracer
#ifdef ALLOW_GENERIC_ADVDIFF
          diagSufx = GAD_DIAG_SUFX( tracerId, myThid )
#else
          diagSufx = 'aaaa'
#endif
          diagName = 'DFrI'//diagSufx
        ELSEIF ( tracerId.EQ. -1 ) THEN
          diagName = 'VISrI_Um'
        ELSEIF ( tracerId.EQ. -2 ) THEN
          diagName = 'VISrI_Vm'
        ELSE
          STOP 'IMPLIDIFF: should never reach this point !'
        ENDIF
        IF ( DIAGNOSTICS_IS_ON(diagName,myThid) ) THEN
         DO k= 1,Nr
          IF ( k.EQ.1 ) THEN
C-  Note: Needs to call DIAGNOSTICS_FILL at level k=1 even if array == 0
C         otherwise counter is not incremented !!
            DO j=1-OLy,sNy+OLy
             DO i=1-OLx,sNx+OLx
               df(i,j) = 0. _d 0
             ENDDO
            ENDDO
          ELSEIF ( tracerId.GE.1 ) THEN
#ifdef TARGET_NEC_SX
            DO j=1-OLy,sNy+OLy
             DO i=1-OLx,sNx+OLx
#else
            DO j=1,sNy
             DO i=1,sNx
#endif
               df(i,j) =
     &             -rA(i,j,bi,bj)*deepFac2F(k)*rhoFacF(k)
     &            * KappaRX(i,j,k)*recip_drC(k)
     &            * (gXnm1(i,j,k,bi,bj) - gXnm1(i,j,k-1,bi,bj))*rkSign
             ENDDO
            ENDDO
          ELSEIF ( tracerId.EQ.-1 ) THEN
#ifdef TARGET_NEC_SX
            DO j=1-OLy,sNy+OLy
             DO i=1-OLx,sNx+OLx
#else
            DO j=1,sNy
             DO i=1,sNx+1
#endif
               df(i,j) =
     &             -rAw(i,j,bi,bj)*deepFac2F(k)*rhoFacF(k)
     &            * KappaRX(i,j,k)*recip_drC(k)
     &            * (gXnm1(i,j,k,bi,bj) - gXnm1(i,j,k-1,bi,bj))*rkSign
     &            * _maskW(i,j,k,bi,bj)
     &            * _maskW(i,j,k-1,bi,bj)
             ENDDO
            ENDDO
          ELSEIF ( tracerId.EQ.-2 ) THEN
#ifdef TARGET_NEC_SX
            DO j=1-OLy,sNy+OLy
             DO i=1-OLx,sNx+OLx
#else
            DO j=1,sNy+1
             DO i=1,sNx
#endif
               df(i,j) =
     &             -rAs(i,j,bi,bj)*deepFac2F(k)*rhoFacF(k)
     &            * KappaRX(i,j,k)*recip_drC(k)
     &            * (gXnm1(i,j,k,bi,bj) - gXnm1(i,j,k-1,bi,bj))*rkSign
     &            * _maskS(i,j,k,bi,bj)
     &            * _maskS(i,j,k-1,bi,bj)
             ENDDO
            ENDDO
          ENDIF
          CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
         ENDDO
        ENDIF
      ENDIF
#endif /* ALLOW_DIAGNOSTICS */

      RETURN
      END
1	mlosch	1.28	C $Header: /u/gcmpack/MITgcm/model/src/impldiff.F,v 1.27 2009/06/26 23:10:09 jahn Exp $
2	cnh	1.17	C $Name: $
3	adcroft	1.1
4	jmc	1.23	#include "PACKAGES_CONFIG.h"
5	cnh	1.7	#include "CPP_OPTIONS.h"
6	adcroft	1.1
7	cnh	1.17	CBOP
8			C !ROUTINE: IMPLDIFF
9			C !INTERFACE:
10	adcroft	1.1	SUBROUTINE IMPLDIFF( bi, bj, iMin, iMax, jMin, jMax,
11	jmc	1.23	I tracerId, KappaRX, recip_hFac,
12	adcroft	1.8	U gXNm1,
13	adcroft	1.1	I myThid )
14	cnh	1.17	C !DESCRIPTION: \bv
15			C ==========================================================
16			C \| S/R IMPLDIFF
17			C \| o Solve implicit diffusion equation for vertical
18			C \| diffusivity.
19			C ==========================================================
20			C \| o Recoded from 2d intermediate fields to 3d to reduce
21			C \| TAMC storage
22			C \| o Fixed missing masks for fields a(), c()
23			C ==========================================================
24			C \ev
25
26			C !USES:
27	cnh	1.5	IMPLICIT NONE
28			C == Global data ==
29	adcroft	1.1	#include "SIZE.h"
30			#include "DYNVARS.h"
31	cnh	1.2	#include "EEPARAMS.h"
32	adcroft	1.1	#include "PARAMS.h"
33			#include "GRID.h"
34	jahn	1.27	#ifdef ALLOW_GENERIC_ADVDIFF
35			#include "GAD.h"
36			#endif
37			#ifdef ALLOW_LONGSTEP
38			#include "LONGSTEP_PARAMS.h"
39			#endif
40			#ifdef ALLOW_PTRACERS
41			#include "PTRACERS_SIZE.h"
42			#include "PTRACERS_PARAMS.h"
43			#endif
44	heimbach	1.9	#ifdef ALLOW_AUTODIFF_TAMC
45			#include "tamc_keys.h"
46			#endif
47
48	cnh	1.17	C !INPUT/OUTPUT PARAMETERS:
49	adcroft	1.1	C == Routine Arguments ==
50	jmc	1.23	C tracerId :: tracer Identificator (if > 0) ;
51			C = 0 or < 0 when solving vertical viscosity implicitly for U or V
52	adcroft	1.1	INTEGER bi,bj,iMin,iMax,jMin,jMax
53	jmc	1.23	INTEGER tracerId
54	adcroft	1.8	_RL KappaRX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
55			_RS recip_hFac(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
56			_RL gXnm1(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
57	adcroft	1.1	INTEGER myThid
58	cnh	1.5
59	cnh	1.17	C !LOCAL VARIABLES:
60	adcroft	1.1	C == Local variables ==
61			INTEGER i,j,k
62	jmc	1.23	_RL deltaTX(Nr)
63	cnh	1.17	_RL gYnm1(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
64	heimbach	1.12	_RL a(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
65			_RL b(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
66			_RL c(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
67			_RL bet(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
68	cnh	1.6	_RL gam(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
69	jmc	1.23	#ifdef ALLOW_DIAGNOSTICS
70			CHARACTER*8 diagName
71			CHARACTER*4 diagSufx
72			#ifdef ALLOW_GENERIC_ADVDIFF
73			CHARACTER*4 GAD_DIAG_SUFX
74			EXTERNAL GAD_DIAG_SUFX
75			#endif
76			LOGICAL DIAGNOSTICS_IS_ON
77			EXTERNAL DIAGNOSTICS_IS_ON
78			_RL df (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
79			#endif /* ALLOW_DIAGNOSTICS */
80	cnh	1.17	CEOP
81	jmc	1.19
82	heimbach	1.21	cph(
83			cph Not good for TAF: may create irreducible control flow graph
84			cph IF (Nr.LE.1) RETURN
85			cph)
86	adcroft	1.1
87	jahn	1.27	#ifdef ALLOW_PTRACERS
88			IF ( tracerId.GE.GAD_TR1) THEN
89			DO k=1,Nr
90			deltaTX(k) = PTRACERS_dTLev(k)
91			ENDDO
92			ELSEIF ( tracerId.GE.1 ) THEN
93			#else
94	jmc	1.23	IF ( tracerId.GE.1 ) THEN
95	jahn	1.27	#endif
96	jmc	1.23	DO k=1,Nr
97			deltaTX(k) = dTtracerLev(k)
98			ENDDO
99			ELSE
100			DO k=1,Nr
101			deltaTX(k) = deltaTmom
102			ENDDO
103			ENDIF
104
105	heimbach	1.12	C-- Initialise
106	heimbach	1.16	DO k=1,Nr
107	mlosch	1.28	#ifdef TARGET_NEC_SX
108			DO j=1-OLy,sNy+OLy
109			DO i=1-OLx,sNx+OLx
110			#else
111	heimbach	1.16	DO j=jMin,jMax
112			DO i=iMin,iMax
113	mlosch	1.28	#endif
114	heimbach	1.16	gYNm1(i,j,k,bi,bj) = 0. _d 0
115			ENDDO
116			ENDDO
117			ENDDO
118	heimbach	1.12
119			C-- Old aLower
120	mlosch	1.28	#ifdef TARGET_NEC_SX
121			DO j=1-OLy,sNy+OLy
122			DO i=1-OLx,sNx+OLx
123			#else
124	heimbach	1.14	DO j=jMin,jMax
125			DO i=iMin,iMax
126	mlosch	1.28	#endif
127	heimbach	1.12	a(i,j,1) = 0. _d 0
128			ENDDO
129			ENDDO
130			DO k=2,Nr
131	mlosch	1.28	#ifdef TARGET_NEC_SX
132			DO j=1-OLy,sNy+OLy
133			DO i=1-OLx,sNx+OLx
134			#else
135	heimbach	1.14	DO j=jMin,jMax
136			DO i=iMin,iMax
137	mlosch	1.28	#endif
138	jmc	1.23	a(i,j,k) = -deltaTX(k)recip_hFac(i,j,k,bi,bj)recip_drF(k)
139	jmc	1.26	& recip_deepFac2C(k)recip_rhoFacC(k)
140	heimbach	1.12	& KappaRX(i,j, k )recip_drC( k )
141	jmc	1.26	& deepFac2F(k)rhoFacF(k)
142	mlosch	1.18	IF (recip_hFac(i,j,k-1,bi,bj).EQ.0.) a(i,j,k)=0.
143	heimbach	1.12	ENDDO
144			ENDDO
145			ENDDO
146
147			C-- Old aUpper
148			DO k=1,Nr-1
149	mlosch	1.28	#ifdef TARGET_NEC_SX
150			DO j=1-OLy,sNy+OLy
151			DO i=1-OLx,sNx+OLx
152			#else
153	heimbach	1.14	DO j=jMin,jMax
154			DO i=iMin,iMax
155	mlosch	1.28	#endif
156	jmc	1.23	c(i,j,k) = -deltaTX(k)recip_hFac(i,j,k,bi,bj)recip_drF(k)
157	jmc	1.26	& recip_deepFac2C(k)recip_rhoFacC(k)
158	heimbach	1.12	& KappaRX(i,j,k+1)recip_drC(k+1)
159	jmc	1.26	& deepFac2F(k+1)rhoFacF(k+1)
160	adcroft	1.13	IF (recip_hFac(i,j,k+1,bi,bj).EQ.0.) c(i,j,k)=0.
161	heimbach	1.12	ENDDO
162			ENDDO
163			ENDDO
164	mlosch	1.28	#ifdef TARGET_NEC_SX
165			DO j=1-OLy,sNy+OLy
166			DO i=1-OLx,sNx+OLx
167			#else
168	heimbach	1.14	DO j=jMin,jMax
169			DO i=iMin,iMax
170	mlosch	1.28	#endif
171	heimbach	1.12	c(i,j,Nr) = 0. _d 0
172			ENDDO
173			ENDDO
174
175			C-- Old aCenter
176			DO k=1,Nr
177	mlosch	1.28	#ifdef TARGET_NEC_SX
178			DO j=1-OLy,sNy+OLy
179			DO i=1-OLx,sNx+OLx
180			#else
181	heimbach	1.14	DO j=jMin,jMax
182			DO i=iMin,iMax
183	mlosch	1.28	#endif
184	heimbach	1.12	b(i,j,k) = 1. _d 0 - c(i,j,k) - a(i,j,k)
185			ENDDO
186			ENDDO
187			ENDDO
188
189			C-- Old and new gam, bet are the same
190			DO k=1,Nr
191	mlosch	1.28	#ifdef TARGET_NEC_SX
192			DO j=1-OLy,sNy+OLy
193			DO i=1-OLx,sNx+OLx
194			#else
195	heimbach	1.14	DO j=jMin,jMax
196			DO i=iMin,iMax
197	mlosch	1.28	#endif
198	jmc	1.22	bet(i,j,k) = 1. _d 0
199	heimbach	1.12	gam(i,j,k) = 0. _d 0
200			ENDDO
201			ENDDO
202			ENDDO
203
204	heimbach	1.10	C-- Only need do anything if Nr>1
205			IF (Nr.GT.1) THEN
206
207	heimbach	1.12	k = 1
208	cnh	1.5	C-- Beginning of forward sweep (top level)
209	mlosch	1.28	#ifdef TARGET_NEC_SX
210			DO j=1-OLy,sNy+OLy
211			DO i=1-OLx,sNx+OLx
212			#else
213	adcroft	1.1	DO j=jMin,jMax
214			DO i=iMin,iMax
215	mlosch	1.28	#endif
216	heimbach	1.12	IF (b(i,j,1).NE.0.) bet(i,j,1) = 1. _d 0 / b(i,j,1)
217	adcroft	1.1	ENDDO
218			ENDDO
219	heimbach	1.10
220	adcroft	1.1	ENDIF
221	heimbach	1.9
222	cnh	1.5	C-- Middle of forward sweep
223	jmc	1.20	IF (Nr.GE.2) THEN
224	heimbach	1.10
225	heimbach	1.12	CADJ loop = sequential
226			DO k=2,Nr
227	heimbach	1.9
228	mlosch	1.28	#ifdef TARGET_NEC_SX
229			DO j=1-OLy,sNy+OLy
230			DO i=1-OLx,sNx+OLx
231			#else
232	adcroft	1.1	DO j=jMin,jMax
233			DO i=iMin,iMax
234	mlosch	1.28	#endif
235	heimbach	1.12	gam(i,j,k) = c(i,j,k-1)*bet(i,j,k-1)
236			IF ( ( b(i,j,k) - a(i,j,k)*gam(i,j,k) ) .NE. 0.)
237			& bet(i,j,k) = 1. _d 0 / ( b(i,j,k) - a(i,j,k)*gam(i,j,k) )
238	adcroft	1.1	ENDDO
239			ENDDO
240	heimbach	1.9
241	adcroft	1.1	ENDDO
242	heimbach	1.10
243	adcroft	1.1	ENDIF
244	heimbach	1.10
245	heimbach	1.11
246	mlosch	1.28	#ifdef TARGET_NEC_SX
247			DO j=1-OLy,sNy+OLy
248			DO i=1-OLx,sNx+OLx
249			#else
250	heimbach	1.12	DO j=jMin,jMax
251			DO i=iMin,iMax
252	mlosch	1.28	#endif
253	heimbach	1.12	gYNm1(i,j,1,bi,bj) = gXNm1(i,j,1,bi,bj)*bet(i,j,1)
254	heimbach	1.10	ENDDO
255	heimbach	1.12	ENDDO
256			DO k=2,Nr
257	mlosch	1.28	#ifdef TARGET_NEC_SX
258			DO j=1-OLy,sNy+OLy
259			DO i=1-OLx,sNx+OLx
260			#else
261	heimbach	1.10	DO j=jMin,jMax
262			DO i=iMin,iMax
263	mlosch	1.28	#endif
264	heimbach	1.12	gYnm1(i,j,k,bi,bj) = bet(i,j,k)*
265			& (gXnm1(i,j,k,bi,bj) - a(i,j,k)*gYnm1(i,j,k-1,bi,bj))
266	heimbach	1.9	ENDDO
267			ENDDO
268	heimbach	1.12	ENDDO
269	heimbach	1.9
270
271	heimbach	1.12	C-- Backward sweep
272			CADJ loop = sequential
273			DO k=Nr-1,1,-1
274	mlosch	1.28	#ifdef TARGET_NEC_SX
275			DO j=1-OLy,sNy+OLy
276			DO i=1-OLx,sNx+OLx
277			#else
278	heimbach	1.12	DO j=jMin,jMax
279			DO i=iMin,iMax
280	mlosch	1.28	#endif
281	heimbach	1.12	gYnm1(i,j,k,bi,bj)=gYnm1(i,j,k,bi,bj)
282			& -gam(i,j,k+1)*gYnm1(i,j,k+1,bi,bj)
283			ENDDO
284	adcroft	1.1	ENDDO
285			ENDDO
286	heimbach	1.9
287	heimbach	1.12	DO k=1,Nr
288	mlosch	1.28	#ifdef TARGET_NEC_SX
289			DO j=1-OLy,sNy+OLy
290			DO i=1-OLx,sNx+OLx
291			#else
292	adcroft	1.1	DO j=jMin,jMax
293			DO i=iMin,iMax
294	mlosch	1.28	#endif
295	heimbach	1.12	gXnm1(i,j,k,bi,bj)=gYnm1(i,j,k,bi,bj)
296	adcroft	1.1	ENDDO
297			ENDDO
298			ENDDO
299
300	jmc	1.23	#ifdef ALLOW_DIAGNOSTICS
301	jmc	1.25	IF ( useDiagnostics .AND.tracerId.NE.0 ) THEN
302			IF ( tracerId.GE. 1 ) THEN
303	jmc	1.23	C-- Set diagnostic suffix for the current tracer
304			#ifdef ALLOW_GENERIC_ADVDIFF
305	jmc	1.25	diagSufx = GAD_DIAG_SUFX( tracerId, myThid )
306	jmc	1.23	#else
307	jmc	1.25	diagSufx = 'aaaa'
308	jmc	1.23	#endif
309	jmc	1.25	diagName = 'DFrI'//diagSufx
310			ELSEIF ( tracerId.EQ. -1 ) THEN
311			diagName = 'VISrI_Um'
312			ELSEIF ( tracerId.EQ. -2 ) THEN
313			diagName = 'VISrI_Vm'
314			ELSE
315			STOP 'IMPLIDIFF: should never reach this point !'
316			ENDIF
317	jmc	1.23	IF ( DIAGNOSTICS_IS_ON(diagName,myThid) ) THEN
318			DO k= 1,Nr
319			IF ( k.EQ.1 ) THEN
320			C- Note: Needs to call DIAGNOSTICS_FILL at level k=1 even if array == 0
321			C otherwise counter is not incremented !!
322			DO j=1-OLy,sNy+OLy
323			DO i=1-OLx,sNx+OLx
324			df(i,j) = 0. _d 0
325			ENDDO
326			ENDDO
327	jmc	1.25	ELSEIF ( tracerId.GE.1 ) THEN
328	mlosch	1.28	#ifdef TARGET_NEC_SX
329			DO j=1-OLy,sNy+OLy
330			DO i=1-OLx,sNx+OLx
331			#else
332	jmc	1.23	DO j=1,sNy
333			DO i=1,sNx
334	mlosch	1.28	#endif
335	jmc	1.23	df(i,j) =
336	jmc	1.26	& -rA(i,j,bi,bj)deepFac2F(k)rhoFacF(k)
337	jmc	1.25	& * KappaRX(i,j,k)*recip_drC(k)
338			& * (gXnm1(i,j,k,bi,bj) - gXnm1(i,j,k-1,bi,bj))*rkSign
339			ENDDO
340			ENDDO
341			ELSEIF ( tracerId.EQ.-1 ) THEN
342	mlosch	1.28	#ifdef TARGET_NEC_SX
343			DO j=1-OLy,sNy+OLy
344			DO i=1-OLx,sNx+OLx
345			#else
346	jmc	1.25	DO j=1,sNy
347			DO i=1,sNx+1
348	mlosch	1.28	#endif
349	jmc	1.25	df(i,j) =
350	jmc	1.26	& -rAw(i,j,bi,bj)deepFac2F(k)rhoFacF(k)
351	jmc	1.25	& * KappaRX(i,j,k)*recip_drC(k)
352			& * (gXnm1(i,j,k,bi,bj) - gXnm1(i,j,k-1,bi,bj))*rkSign
353			& * _maskW(i,j,k,bi,bj)
354			& * _maskW(i,j,k-1,bi,bj)
355			ENDDO
356			ENDDO
357			ELSEIF ( tracerId.EQ.-2 ) THEN
358	mlosch	1.28	#ifdef TARGET_NEC_SX
359			DO j=1-OLy,sNy+OLy
360			DO i=1-OLx,sNx+OLx
361			#else
362	jmc	1.25	DO j=1,sNy+1
363			DO i=1,sNx
364	mlosch	1.28	#endif
365	jmc	1.25	df(i,j) =
366	jmc	1.26	& -rAs(i,j,bi,bj)deepFac2F(k)rhoFacF(k)
367	jmc	1.23	& * KappaRX(i,j,k)*recip_drC(k)
368	jmc	1.25	& * (gXnm1(i,j,k,bi,bj) - gXnm1(i,j,k-1,bi,bj))*rkSign
369			& * _maskS(i,j,k,bi,bj)
370			& * _maskS(i,j,k-1,bi,bj)
371	jmc	1.23	ENDDO
372			ENDDO
373			ENDIF
374	jmc	1.24	CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
375	jmc	1.23	ENDDO
376			ENDIF
377			ENDIF
378			#endif /* ALLOW_DIAGNOSTICS */
379
380	adcroft	1.1	RETURN
381			END