model/src/impldiff.F

C $Header: /u/gcmpack/models/MITgcmUV/model/src/impldiff.F,v 1.8 1999/05/18 18:01:13 adcroft 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_impl3d, key = idkey
CADJ store gXNm1 = comlev1_impl2d, 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 = comlev1_impl3d, key = kkey
CADJ store bet  = comlev1_impl3d, 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_impl3d, key = kkey
CADJ store gXNm1(:,:,k-1:k,bi,bj) = comlev1_impl3d, 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_impl2d, key = idkey
CADJ store a,b  = comlev1_impl2d, key = idkey
CADJ store bet  = comlev1_impl2d, 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 gXnm1 = comlev1_impl2d, 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_impl2d, 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.9	C $Header: /u/gcmpack/models/MITgcmUV/model/src/impldiff.F,v 1.8 1999/05/18 18:01:13 adcroft 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			CADJ store bet = comlev1_impl3d, key = idkey
65			CADJ store gXNm1 = comlev1_impl2d, key = idkey
66			#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			CADJ store ckm1 = comlev1_impl3d, key = kkey
96			CADJ store bet = comlev1_impl3d, key = kkey
97			#endif
98
99	adcroft	1.1	DO j=jMin,jMax
100			DO i=iMin,iMax
101			gam(i,j,k)=ckm1(i,j)*bet(i,j)
102	heimbach	1.10	ENDDO
103			ENDDO
104
105			DO j=jMin,jMax
106			DO i=iMin,iMax
107	adcroft	1.1	bet(i,j)=b(i,j)-a(i,j)*gam(i,j,k)
108			IF (bet(i,j).NE.0.) bet(i,j)=1. / bet(i,j)
109	heimbach	1.9	ENDDO
110			ENDDO
111
112			#ifdef ALLOW_AUTODIFF_TAMC
113			CADJ store bet = comlev1_impl3d, key = kkey
114			CADJ store gXNm1(:,:,k-1:k,bi,bj) = comlev1_impl3d, key = kkey
115			#endif
116
117			DO j=jMin,jMax
118			DO i=iMin,iMax
119	adcroft	1.8	gXnm1(i,j,k,bi,bj)=(gXnm1(i,j,k,bi,bj)
120			& -a(i,j)gXnm1(i,j,k-1,bi,bj))bet(i,j)
121	adcroft	1.1	ENDDO
122			ENDDO
123	heimbach	1.9
124	adcroft	1.1	ENDDO
125	heimbach	1.10
126	adcroft	1.1	ENDIF
127	heimbach	1.10
128	cnh	1.6	IF (Nr.GT.1) THEN
129	cnh	1.5	C-- End of forward sweep (bottom level)
130	adcroft	1.1	DO j=jMin,jMax
131			DO i=iMin,iMax
132			ckm1(i,j)=c(i,j)
133	adcroft	1.8	a(i,j)=-deltaTXrecip_hFac(i,j,Nr,bi,bj)recip_drF(Nr)
134			& KappaRX(i,j, Nr )recip_drC( Nr )
135	adcroft	1.1	b(i,j)=1.-a(i,j)
136			ENDDO
137			ENDDO
138	heimbach	1.9
139			#ifdef ALLOW_AUTODIFF_TAMC
140			CADJ store ckm1 = comlev1_impl2d, key = idkey
141			CADJ store a,b = comlev1_impl2d, key = idkey
142			CADJ store bet = comlev1_impl2d, key = idkey
143			#endif
144
145	adcroft	1.1	DO j=jMin,jMax
146			DO i=iMin,iMax
147	cnh	1.6	gam(i,j,Nr)=ckm1(i,j)*bet(i,j)
148	heimbach	1.10	ENDDO
149			ENDDO
150			DO j=jMin,jMax
151			DO i=iMin,iMax
152	cnh	1.6	bet(i,j)=b(i,j)-a(i,j)*gam(i,j,Nr)
153	adcroft	1.1	IF (bet(i,j).NE.0.) bet(i,j)=1. / bet(i,j)
154	heimbach	1.9	ENDDO
155			ENDDO
156
157			#ifdef ALLOW_AUTODIFF_TAMC
158			CADJ store gXnm1 = comlev1_impl2d, key = idkey
159			#endif
160
161			DO j=jMin,jMax
162			DO i=iMin,iMax
163	adcroft	1.8	gXnm1(i,j,Nr,bi,bj)=(gXnm1(i,j,Nr,bi,bj)
164			& -a(i,j)gXnm1(i,j,Nr-1,bi,bj))bet(i,j)
165	adcroft	1.1	ENDDO
166			ENDDO
167	heimbach	1.9
168			#ifdef ALLOW_AUTODIFF_TAMC
169			CADJ store gam = comlev1_impl2d, key = idkey
170			#endif
171
172	cnh	1.5	C-- Backward sweep
173	cnh	1.6	DO k=Nr-1,1,-1
174	adcroft	1.1	DO j=jMin,jMax
175			DO i=iMin,iMax
176	adcroft	1.8	gXnm1(i,j,k,bi,bj)=gXnm1(i,j,k,bi,bj)
177			& -gam(i,j,k+1)*gXnm1(i,j,k+1,bi,bj)
178	adcroft	1.1	ENDDO
179			ENDDO
180			ENDDO
181	heimbach	1.10
182	adcroft	1.1	ENDIF
183
184			RETURN
185			END