model/src/impldiff.F

C $Header: /u/gcmpack/MITgcm/model/src/impldiff.F,v 1.22 2004/10/22 20:14:46 jmc 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_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
      INTEGER     kk
      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)

      IF ( tracerId.GE.1 ) THEN
        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
       DO j=jMin,jMax
        DO i=iMin,iMax
         gYNm1(i,j,k,bi,bj) = 0. _d 0
        ENDDO
       ENDDO
      ENDDO

C--   Old aLower
      DO j=jMin,jMax
       DO i=iMin,iMax
         a(i,j,1) = 0. _d 0 
       ENDDO
      ENDDO
      DO k=2,Nr
       DO j=jMin,jMax
        DO i=iMin,iMax
          a(i,j,k) = -deltaTX(k)*recip_hFac(i,j,k,bi,bj)*recip_drF(k)
     &               *KappaRX(i,j, k )*recip_drC( 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
       DO j=jMin,jMax
        DO i=iMin,iMax
          c(i,j,k) = -deltaTX(k)*recip_hFac(i,j,k,bi,bj)*recip_drF(k)
     &               *KappaRX(i,j,k+1)*recip_drC(k+1)
          IF (recip_hFac(i,j,k+1,bi,bj).EQ.0.) c(i,j,k)=0.
        ENDDO
       ENDDO
      ENDDO
      DO j=jMin,jMax
       DO i=iMin,iMax
         c(i,j,Nr) = 0. _d 0
       ENDDO
      ENDDO

C--   Old aCenter
      DO k=1,Nr
       DO j=jMin,jMax
        DO i=iMin,iMax
          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
       DO j=jMin,jMax
        DO i=iMin,iMax
          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)
       DO j=jMin,jMax
        DO i=iMin,iMax
         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

        DO j=jMin,jMax
         DO i=iMin,iMax
          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


      DO j=jMin,jMax
       DO i=iMin,iMax
        gYNm1(i,j,1,bi,bj) = gXNm1(i,j,1,bi,bj)*bet(i,j,1)
       ENDDO
      ENDDO
      DO k=2,Nr
       DO j=jMin,jMax
        DO i=iMin,iMax
         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
        DO j=jMin,jMax
         DO i=iMin,iMax
          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
        DO j=jMin,jMax
         DO i=iMin,iMax
          gXnm1(i,j,k,bi,bj)=gYnm1(i,j,k,bi,bj)
         ENDDO
        ENDDO
       ENDDO

#ifdef ALLOW_DIAGNOSTICS
      IF ( useDiagnostics .AND.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
        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
          ELSE
            DO j=1,sNy
             DO i=1,sNx
               df(i,j) =
     &              rA(i,j,bi,bj)
     &            * KappaRX(i,j,k)*recip_drC(k)
     &            * (gXnm1(i,j,k,bi,bj) - gXnm1(i,j,k-1,bi,bj))
             ENDDO
            ENDDO
          ENDIF
          kk = -k
          CALL DIAGNOSTICS_FILL(df,diagName, kk,1, 2,bi,bj, myThid)
         ENDDO
        ENDIF
      ENDIF
#endif /* ALLOW_DIAGNOSTICS */

      RETURN
      END
1	jmc	1.23	C $Header: /u/gcmpack/MITgcm/model/src/impldiff.F,v 1.22 2004/10/22 20:14:46 jmc 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	heimbach	1.9	#ifdef ALLOW_AUTODIFF_TAMC
35			#include "tamc_keys.h"
36			#endif
37
38	cnh	1.17	C !INPUT/OUTPUT PARAMETERS:
39	adcroft	1.1	C == Routine Arguments ==
40	jmc	1.23	C tracerId :: tracer Identificator (if > 0) ;
41			C = 0 or < 0 when solving vertical viscosity implicitly for U or V
42	adcroft	1.1	INTEGER bi,bj,iMin,iMax,jMin,jMax
43	jmc	1.23	INTEGER tracerId
44	adcroft	1.8	_RL KappaRX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
45			_RS recip_hFac(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
46			_RL gXnm1(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
47	adcroft	1.1	INTEGER myThid
48	cnh	1.5
49	cnh	1.17	C !LOCAL VARIABLES:
50	adcroft	1.1	C == Local variables ==
51			INTEGER i,j,k
52	jmc	1.23	_RL deltaTX(Nr)
53	cnh	1.17	_RL gYnm1(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
54	heimbach	1.12	_RL a(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
55			_RL b(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
56			_RL c(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
57			_RL bet(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
58	cnh	1.6	_RL gam(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
59	jmc	1.23	#ifdef ALLOW_DIAGNOSTICS
60			INTEGER kk
61			CHARACTER*8 diagName
62			CHARACTER*4 diagSufx
63			#ifdef ALLOW_GENERIC_ADVDIFF
64			CHARACTER*4 GAD_DIAG_SUFX
65			EXTERNAL GAD_DIAG_SUFX
66			#endif
67			LOGICAL DIAGNOSTICS_IS_ON
68			EXTERNAL DIAGNOSTICS_IS_ON
69			_RL df (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
70			#endif /* ALLOW_DIAGNOSTICS */
71	cnh	1.17	CEOP
72	jmc	1.19
73	heimbach	1.21	cph(
74			cph Not good for TAF: may create irreducible control flow graph
75			cph IF (Nr.LE.1) RETURN
76			cph)
77	adcroft	1.1
78	jmc	1.23	IF ( tracerId.GE.1 ) THEN
79			DO k=1,Nr
80			deltaTX(k) = dTtracerLev(k)
81			ENDDO
82			ELSE
83			DO k=1,Nr
84			deltaTX(k) = deltaTmom
85			ENDDO
86			ENDIF
87
88	heimbach	1.12	C-- Initialise
89	heimbach	1.16	DO k=1,Nr
90			DO j=jMin,jMax
91			DO i=iMin,iMax
92			gYNm1(i,j,k,bi,bj) = 0. _d 0
93			ENDDO
94			ENDDO
95			ENDDO
96	heimbach	1.12
97			C-- Old aLower
98	heimbach	1.14	DO j=jMin,jMax
99			DO i=iMin,iMax
100	heimbach	1.12	a(i,j,1) = 0. _d 0
101			ENDDO
102			ENDDO
103			DO k=2,Nr
104	heimbach	1.14	DO j=jMin,jMax
105			DO i=iMin,iMax
106	jmc	1.23	a(i,j,k) = -deltaTX(k)recip_hFac(i,j,k,bi,bj)recip_drF(k)
107	heimbach	1.12	& KappaRX(i,j, k )recip_drC( k )
108	mlosch	1.18	IF (recip_hFac(i,j,k-1,bi,bj).EQ.0.) a(i,j,k)=0.
109	heimbach	1.12	ENDDO
110			ENDDO
111			ENDDO
112
113			C-- Old aUpper
114			DO k=1,Nr-1
115	heimbach	1.14	DO j=jMin,jMax
116			DO i=iMin,iMax
117	jmc	1.23	c(i,j,k) = -deltaTX(k)recip_hFac(i,j,k,bi,bj)recip_drF(k)
118	heimbach	1.12	& KappaRX(i,j,k+1)recip_drC(k+1)
119	adcroft	1.13	IF (recip_hFac(i,j,k+1,bi,bj).EQ.0.) c(i,j,k)=0.
120	heimbach	1.12	ENDDO
121			ENDDO
122			ENDDO
123	heimbach	1.14	DO j=jMin,jMax
124			DO i=iMin,iMax
125	heimbach	1.12	c(i,j,Nr) = 0. _d 0
126			ENDDO
127			ENDDO
128
129			C-- Old aCenter
130			DO k=1,Nr
131	heimbach	1.14	DO j=jMin,jMax
132			DO i=iMin,iMax
133	heimbach	1.12	b(i,j,k) = 1. _d 0 - c(i,j,k) - a(i,j,k)
134			ENDDO
135			ENDDO
136			ENDDO
137
138			C-- Old and new gam, bet are the same
139			DO k=1,Nr
140	heimbach	1.14	DO j=jMin,jMax
141			DO i=iMin,iMax
142	jmc	1.22	bet(i,j,k) = 1. _d 0
143	heimbach	1.12	gam(i,j,k) = 0. _d 0
144			ENDDO
145			ENDDO
146			ENDDO
147
148	heimbach	1.10	C-- Only need do anything if Nr>1
149			IF (Nr.GT.1) THEN
150
151	heimbach	1.12	k = 1
152	cnh	1.5	C-- Beginning of forward sweep (top level)
153	adcroft	1.1	DO j=jMin,jMax
154			DO i=iMin,iMax
155	heimbach	1.12	IF (b(i,j,1).NE.0.) bet(i,j,1) = 1. _d 0 / b(i,j,1)
156	adcroft	1.1	ENDDO
157			ENDDO
158	heimbach	1.10
159	adcroft	1.1	ENDIF
160	heimbach	1.9
161	cnh	1.5	C-- Middle of forward sweep
162	jmc	1.20	IF (Nr.GE.2) THEN
163	heimbach	1.10
164	heimbach	1.12	CADJ loop = sequential
165			DO k=2,Nr
166	heimbach	1.9
167	adcroft	1.1	DO j=jMin,jMax
168			DO i=iMin,iMax
169	heimbach	1.12	gam(i,j,k) = c(i,j,k-1)*bet(i,j,k-1)
170			IF ( ( b(i,j,k) - a(i,j,k)*gam(i,j,k) ) .NE. 0.)
171			& bet(i,j,k) = 1. _d 0 / ( b(i,j,k) - a(i,j,k)*gam(i,j,k) )
172	adcroft	1.1	ENDDO
173			ENDDO
174	heimbach	1.9
175	adcroft	1.1	ENDDO
176	heimbach	1.10
177	adcroft	1.1	ENDIF
178	heimbach	1.10
179	heimbach	1.11
180	heimbach	1.12	DO j=jMin,jMax
181			DO i=iMin,iMax
182			gYNm1(i,j,1,bi,bj) = gXNm1(i,j,1,bi,bj)*bet(i,j,1)
183	heimbach	1.10	ENDDO
184	heimbach	1.12	ENDDO
185			DO k=2,Nr
186	heimbach	1.10	DO j=jMin,jMax
187			DO i=iMin,iMax
188	heimbach	1.12	gYnm1(i,j,k,bi,bj) = bet(i,j,k)*
189			& (gXnm1(i,j,k,bi,bj) - a(i,j,k)*gYnm1(i,j,k-1,bi,bj))
190	heimbach	1.9	ENDDO
191			ENDDO
192	heimbach	1.12	ENDDO
193	heimbach	1.9
194
195	heimbach	1.12	C-- Backward sweep
196			CADJ loop = sequential
197			DO k=Nr-1,1,-1
198			DO j=jMin,jMax
199			DO i=iMin,iMax
200			gYnm1(i,j,k,bi,bj)=gYnm1(i,j,k,bi,bj)
201			& -gam(i,j,k+1)*gYnm1(i,j,k+1,bi,bj)
202			ENDDO
203	adcroft	1.1	ENDDO
204			ENDDO
205	heimbach	1.9
206	heimbach	1.12	DO k=1,Nr
207	adcroft	1.1	DO j=jMin,jMax
208			DO i=iMin,iMax
209	heimbach	1.12	gXnm1(i,j,k,bi,bj)=gYnm1(i,j,k,bi,bj)
210	adcroft	1.1	ENDDO
211			ENDDO
212			ENDDO
213
214	jmc	1.23	#ifdef ALLOW_DIAGNOSTICS
215			IF ( useDiagnostics .AND.tracerId.GE.1 ) THEN
216			C-- Set diagnostic suffix for the current tracer
217			#ifdef ALLOW_GENERIC_ADVDIFF
218			diagSufx = GAD_DIAG_SUFX( tracerId, myThid )
219			#else
220			diagSufx = 'aaaa'
221			#endif
222			diagName = 'DFrI'//diagSufx
223			IF ( DIAGNOSTICS_IS_ON(diagName,myThid) ) THEN
224			DO k= 1,Nr
225			IF ( k.EQ.1 ) THEN
226			C- Note: Needs to call DIAGNOSTICS_FILL at level k=1 even if array == 0
227			C otherwise counter is not incremented !!
228			DO j=1-OLy,sNy+OLy
229			DO i=1-OLx,sNx+OLx
230			df(i,j) = 0. _d 0
231			ENDDO
232			ENDDO
233			ELSE
234			DO j=1,sNy
235			DO i=1,sNx
236			df(i,j) =
237			& rA(i,j,bi,bj)
238			& * KappaRX(i,j,k)*recip_drC(k)
239			& * (gXnm1(i,j,k,bi,bj) - gXnm1(i,j,k-1,bi,bj))
240			ENDDO
241			ENDDO
242			ENDIF
243			kk = -k
244			CALL DIAGNOSTICS_FILL(df,diagName, kk,1, 2,bi,bj, myThid)
245			ENDDO
246			ENDIF
247			ENDIF
248			#endif /* ALLOW_DIAGNOSTICS */
249
250	adcroft	1.1	RETURN
251			END