model/src/impldiff.F

C $Header: /u/gcmpack/MITgcm/model/src/impldiff.F,v 1.23 2004/12/16 23:20:06 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
      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
          CALL DIAGNOSTICS_FILL(df,diagName, k,1, 2,bi,bj, myThid)
         ENDDO
        ENDIF
      ENDIF
#endif /* ALLOW_DIAGNOSTICS */

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