model/src/impldiff.F

C $Header: /u/gcmpack/models/MITgcmUV/model/src/impldiff.F,v 1.10 2000/06/12 13:38:50 heimbach Exp $

#include "CPP_OPTIONS.h"

C     /==========================================================\
C     | S/R IMPLDIFF                                             |
C     | o Solve implicit diffusion equation for vertical         |
C     |   diffusivity.                                           |
C     \==========================================================/
      SUBROUTINE IMPLDIFF( bi, bj, iMin, iMax, jMin, jMax,
     I                     deltaTX,KappaRX,recip_hFac,
     U                     gXNm1,
     I                     myThid )
      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     == Routine Arguments ==
      INTEGER bi,bj,iMin,iMax,jMin,jMax
      _RL deltaTX
      _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 ==
      INTEGER i,j,k
      _RL a(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      _RL b(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      _RL c(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      _RL ckm1(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      _RL bet(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      _RL gam(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)

#ifdef ALLOW_AUTODIFF_TAMC
      INTEGER kkey
#endif

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

C--    Beginning of forward sweep (top level)
       DO j=jMin,jMax
        DO i=iMin,iMax
         c(i,j)=-deltaTX*recip_hFac(i,j,1,bi,bj)*recip_drF(1)
     &         *KappaRX(i,j,2)*recip_drC(2)
         b(i,j)=1.-c(i,j)
         bet(i,j)=0.
         IF (b(i,j).NE.0.) bet(i,j)=1. / b(i,j)
        ENDDO
       ENDDO

      ENDIF

#ifdef ALLOW_AUTODIFF_TAMC
CADJ store bet                 = comlev1_impl, key = idkey
CADJ store gXNm1(:,:,:,bi,bj)  = comlev1_impl, key = idkey
#endif

       DO j=jMin,jMax
        DO i=iMin,iMax
         gXNm1(i,j,1,bi,bj) = gXNm1(i,j,1,bi,bj)*bet(i,j)
        ENDDO
       ENDDO

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

       DO k=2,Nr-1

#ifdef ALLOW_AUTODIFF_TAMC
        kkey = (idkey-1)*(Nr-2) + k-1
#endif

        DO j=jMin,jMax
         DO i=iMin,iMax
          ckm1(i,j)=c(i,j)
          a(i,j)=-deltaTX*recip_hFac(i,j,k,bi,bj)*recip_drF(k)
     &         *KappaRX(i,j, k )*recip_drC( k )
          c(i,j)=-deltaTX*recip_hFac(i,j,k,bi,bj)*recip_drF(k)
     &         *KappaRX(i,j,k+1)*recip_drC(k+1)
          b(i,j)=1.-c(i,j)-a(i,j)
         ENDDO
        ENDDO

#ifdef ALLOW_AUTODIFF_TAMC
CADJ store ckm1, bet = comlev1_impl_k, key = kkey
#endif

        DO j=jMin,jMax
         DO i=iMin,iMax
          gam(i,j,k)=ckm1(i,j)*bet(i,j)
         ENDDO
        ENDDO

        DO j=jMin,jMax
         DO i=iMin,iMax
          bet(i,j)=b(i,j)-a(i,j)*gam(i,j,k)
          IF (bet(i,j).NE.0.) bet(i,j)=1. / bet(i,j)
         ENDDO
        ENDDO

#ifdef ALLOW_AUTODIFF_TAMC
CADJ store bet                    = comlev1_impl_k, key = kkey
CADJ store gXNm1(:,:,k-1:k,bi,bj) = comlev1_impl_k, key = kkey
#endif

        DO j=jMin,jMax
         DO i=iMin,iMax
          gXnm1(i,j,k,bi,bj)=(gXnm1(i,j,k,bi,bj)
     &       -a(i,j)*gXnm1(i,j,k-1,bi,bj))*bet(i,j)
         ENDDO
        ENDDO

       ENDDO

      ENDIF


      IF (Nr.GT.1) THEN

C--    End of forward sweep (bottom level)
       DO j=jMin,jMax
        DO i=iMin,iMax
         ckm1(i,j)=c(i,j)
         a(i,j)=-deltaTX*recip_hFac(i,j,Nr,bi,bj)*recip_drF(Nr)
     &        *KappaRX(i,j, Nr )*recip_drC( Nr )
         b(i,j)=1.-a(i,j)
        ENDDO
       ENDDO

#ifdef ALLOW_AUTODIFF_TAMC
CADJ store ckm1 = comlev1_impl, key = idkey
CADJ store a,b  = comlev1_impl, key = idkey
CADJ store bet  = comlev1_impl, key = idkey
#endif

       DO j=jMin,jMax
        DO i=iMin,iMax
         gam(i,j,Nr)=ckm1(i,j)*bet(i,j)
        ENDDO
       ENDDO
       DO j=jMin,jMax
        DO i=iMin,iMax
         bet(i,j)=b(i,j)-a(i,j)*gam(i,j,Nr)
         IF (bet(i,j).NE.0.) bet(i,j)=1. / bet(i,j)
        ENDDO
       ENDDO

#ifdef ALLOW_AUTODIFF_TAMC
CADJ store a,bet              = comlev1_impl, key = idkey
CADJ store gXnm1(:,:,:,bi,bj) = comlev1_impl, key = idkey
#endif

       DO j=jMin,jMax
        DO i=iMin,iMax
         gXnm1(i,j,Nr,bi,bj)=(gXnm1(i,j,Nr,bi,bj)
     &      -a(i,j)*gXnm1(i,j,Nr-1,bi,bj))*bet(i,j)
        ENDDO
       ENDDO

#ifdef ALLOW_AUTODIFF_TAMC
CADJ store gam   = comlev1_impl, key = idkey
#endif

C--    Backward sweep
       DO k=Nr-1,1,-1
        DO j=jMin,jMax
         DO i=iMin,iMax
          gXnm1(i,j,k,bi,bj)=gXnm1(i,j,k,bi,bj)
     &              -gam(i,j,k+1)*gXnm1(i,j,k+1,bi,bj)
         ENDDO
        ENDDO
       ENDDO

      ENDIF

      RETURN
      END
1	heimbach	1.11	C $Header: /u/gcmpack/models/MITgcmUV/model/src/impldiff.F,v 1.10 2000/06/12 13:38:50 heimbach Exp $
2	adcroft	1.1
3	cnh	1.7	#include "CPP_OPTIONS.h"
4	adcroft	1.1
5			C /==========================================================\
6			C \| S/R IMPLDIFF \|
7	cnh	1.5	C \| o Solve implicit diffusion equation for vertical \|
8			C \| diffusivity. \|
9	adcroft	1.1	C \==========================================================/
10			SUBROUTINE IMPLDIFF( bi, bj, iMin, iMax, jMin, jMax,
11	adcroft	1.8	I deltaTX,KappaRX,recip_hFac,
12			U gXNm1,
13	adcroft	1.1	I myThid )
14	cnh	1.5	IMPLICIT NONE
15			C == Global data ==
16	adcroft	1.1	#include "SIZE.h"
17			#include "DYNVARS.h"
18	cnh	1.2	#include "EEPARAMS.h"
19	adcroft	1.1	#include "PARAMS.h"
20			#include "GRID.h"
21	cnh	1.5
22	heimbach	1.9	#ifdef ALLOW_AUTODIFF_TAMC
23			#include "tamc_keys.h"
24			#endif
25
26	adcroft	1.1	C == Routine Arguments ==
27			INTEGER bi,bj,iMin,iMax,jMin,jMax
28	adcroft	1.8	_RL deltaTX
29			_RL KappaRX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
30			_RS recip_hFac(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
31			_RL gXnm1(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
32	adcroft	1.1	INTEGER myThid
33	cnh	1.5
34	adcroft	1.1	C == Local variables ==
35			INTEGER i,j,k
36			_RL a(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
37			_RL b(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
38			_RL c(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
39			_RL ckm1(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
40			_RL bet(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
41	cnh	1.6	_RL gam(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
42	adcroft	1.1
43	heimbach	1.9	#ifdef ALLOW_AUTODIFF_TAMC
44			INTEGER kkey
45			#endif
46
47	heimbach	1.10	C-- Only need do anything if Nr>1
48			IF (Nr.GT.1) THEN
49
50	cnh	1.5	C-- Beginning of forward sweep (top level)
51	adcroft	1.1	DO j=jMin,jMax
52			DO i=iMin,iMax
53	adcroft	1.8	c(i,j)=-deltaTXrecip_hFac(i,j,1,bi,bj)recip_drF(1)
54			& KappaRX(i,j,2)recip_drC(2)
55	adcroft	1.1	b(i,j)=1.-c(i,j)
56			bet(i,j)=0.
57			IF (b(i,j).NE.0.) bet(i,j)=1. / b(i,j)
58			ENDDO
59			ENDDO
60	heimbach	1.10
61	adcroft	1.1	ENDIF
62	heimbach	1.9
63			#ifdef ALLOW_AUTODIFF_TAMC
64	heimbach	1.11	CADJ store bet = comlev1_impl, key = idkey
65			CADJ store gXNm1(:,:,:,bi,bj) = comlev1_impl, key = idkey
66	heimbach	1.9	#endif
67
68			DO j=jMin,jMax
69			DO i=iMin,iMax
70			gXNm1(i,j,1,bi,bj) = gXNm1(i,j,1,bi,bj)*bet(i,j)
71			ENDDO
72			ENDDO
73
74	cnh	1.5	C-- Middle of forward sweep
75	cnh	1.6	IF (Nr.GT.2) THEN
76	heimbach	1.10
77	cnh	1.6	DO k=2,Nr-1
78	heimbach	1.9
79			#ifdef ALLOW_AUTODIFF_TAMC
80			kkey = (idkey-1)*(Nr-2) + k-1
81			#endif
82
83	adcroft	1.1	DO j=jMin,jMax
84			DO i=iMin,iMax
85			ckm1(i,j)=c(i,j)
86	adcroft	1.8	a(i,j)=-deltaTXrecip_hFac(i,j,k,bi,bj)recip_drF(k)
87			& KappaRX(i,j, k )recip_drC( k )
88			c(i,j)=-deltaTXrecip_hFac(i,j,k,bi,bj)recip_drF(k)
89			& KappaRX(i,j,k+1)recip_drC(k+1)
90	adcroft	1.1	b(i,j)=1.-c(i,j)-a(i,j)
91			ENDDO
92			ENDDO
93	heimbach	1.9
94			#ifdef ALLOW_AUTODIFF_TAMC
95	heimbach	1.11	CADJ store ckm1, bet = comlev1_impl_k, key = kkey
96	heimbach	1.9	#endif
97
98	adcroft	1.1	DO j=jMin,jMax
99			DO i=iMin,iMax
100			gam(i,j,k)=ckm1(i,j)*bet(i,j)
101	heimbach	1.10	ENDDO
102			ENDDO
103
104			DO j=jMin,jMax
105			DO i=iMin,iMax
106	adcroft	1.1	bet(i,j)=b(i,j)-a(i,j)*gam(i,j,k)
107			IF (bet(i,j).NE.0.) bet(i,j)=1. / bet(i,j)
108	heimbach	1.9	ENDDO
109			ENDDO
110
111			#ifdef ALLOW_AUTODIFF_TAMC
112	heimbach	1.11	CADJ store bet = comlev1_impl_k, key = kkey
113			CADJ store gXNm1(:,:,k-1:k,bi,bj) = comlev1_impl_k, key = kkey
114	heimbach	1.9	#endif
115
116			DO j=jMin,jMax
117			DO i=iMin,iMax
118	adcroft	1.8	gXnm1(i,j,k,bi,bj)=(gXnm1(i,j,k,bi,bj)
119			& -a(i,j)gXnm1(i,j,k-1,bi,bj))bet(i,j)
120	adcroft	1.1	ENDDO
121			ENDDO
122	heimbach	1.9
123	adcroft	1.1	ENDDO
124	heimbach	1.10
125	adcroft	1.1	ENDIF
126	heimbach	1.10
127	heimbach	1.11
128	cnh	1.6	IF (Nr.GT.1) THEN
129	heimbach	1.11
130	cnh	1.5	C-- End of forward sweep (bottom level)
131	adcroft	1.1	DO j=jMin,jMax
132			DO i=iMin,iMax
133			ckm1(i,j)=c(i,j)
134	adcroft	1.8	a(i,j)=-deltaTXrecip_hFac(i,j,Nr,bi,bj)recip_drF(Nr)
135			& KappaRX(i,j, Nr )recip_drC( Nr )
136	adcroft	1.1	b(i,j)=1.-a(i,j)
137			ENDDO
138			ENDDO
139	heimbach	1.9
140			#ifdef ALLOW_AUTODIFF_TAMC
141	heimbach	1.11	CADJ store ckm1 = comlev1_impl, key = idkey
142			CADJ store a,b = comlev1_impl, key = idkey
143			CADJ store bet = comlev1_impl, key = idkey
144	heimbach	1.9	#endif
145
146	adcroft	1.1	DO j=jMin,jMax
147			DO i=iMin,iMax
148	cnh	1.6	gam(i,j,Nr)=ckm1(i,j)*bet(i,j)
149	heimbach	1.10	ENDDO
150			ENDDO
151			DO j=jMin,jMax
152			DO i=iMin,iMax
153	cnh	1.6	bet(i,j)=b(i,j)-a(i,j)*gam(i,j,Nr)
154	adcroft	1.1	IF (bet(i,j).NE.0.) bet(i,j)=1. / bet(i,j)
155	heimbach	1.9	ENDDO
156			ENDDO
157
158			#ifdef ALLOW_AUTODIFF_TAMC
159	heimbach	1.11	CADJ store a,bet = comlev1_impl, key = idkey
160			CADJ store gXnm1(:,:,:,bi,bj) = comlev1_impl, key = idkey
161	heimbach	1.9	#endif
162
163			DO j=jMin,jMax
164			DO i=iMin,iMax
165	adcroft	1.8	gXnm1(i,j,Nr,bi,bj)=(gXnm1(i,j,Nr,bi,bj)
166			& -a(i,j)gXnm1(i,j,Nr-1,bi,bj))bet(i,j)
167	adcroft	1.1	ENDDO
168			ENDDO
169	heimbach	1.9
170			#ifdef ALLOW_AUTODIFF_TAMC
171	heimbach	1.11	CADJ store gam = comlev1_impl, key = idkey
172	heimbach	1.9	#endif
173
174	cnh	1.5	C-- Backward sweep
175	cnh	1.6	DO k=Nr-1,1,-1
176	adcroft	1.1	DO j=jMin,jMax
177			DO i=iMin,iMax
178	adcroft	1.8	gXnm1(i,j,k,bi,bj)=gXnm1(i,j,k,bi,bj)
179			& -gam(i,j,k+1)*gXnm1(i,j,k+1,bi,bj)
180	adcroft	1.1	ENDDO
181			ENDDO
182			ENDDO
183	heimbach	1.10
184	adcroft	1.1	ENDIF
185
186			RETURN
187			END