model/src/impldiff.F

C $Header: /u/gcmpack/models/MITgcmUV/model/src/impldiff.F,v 1.7 1998/11/06 22:44:46 cnh 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"

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)

      IF (Nr.GT.1) THEN         ! Only need do anything if Nr>1
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)
         gXnm1(i,j,1,bi,bj)=gXnm1(i,j,1,bi,bj)*bet(i,j)
        ENDDO
       ENDDO
      ENDIF
C--   Middle of forward sweep
      IF (Nr.GT.2) THEN
       DO k=2,Nr-1
        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
        DO j=jMin,jMax
         DO i=iMin,iMax
          gam(i,j,k)=ckm1(i,j)*bet(i,j)
          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)
          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
       DO j=jMin,jMax
        DO i=iMin,iMax
         gam(i,j,Nr)=ckm1(i,j)*bet(i,j)
         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)
         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
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	adcroft	1.8	C $Header: /u/gcmpack/models/MITgcmUV/model/src/impldiff.F,v 1.7 1998/11/06 22:44:46 cnh 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	adcroft	1.1	C == Routine Arguments ==
23			INTEGER bi,bj,iMin,iMax,jMin,jMax
24	adcroft	1.8	_RL deltaTX
25			_RL KappaRX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
26			_RS recip_hFac(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
27			_RL gXnm1(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
28	adcroft	1.1	INTEGER myThid
29	cnh	1.5
30	adcroft	1.1	C == Local variables ==
31			INTEGER i,j,k
32			_RL a(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
33			_RL b(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
34			_RL c(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
35			_RL ckm1(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
36			_RL bet(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
37	cnh	1.6	_RL gam(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
38	adcroft	1.1
39	cnh	1.6	IF (Nr.GT.1) THEN ! Only need do anything if Nr>1
40	cnh	1.5	C-- Beginning of forward sweep (top level)
41	adcroft	1.1	DO j=jMin,jMax
42			DO i=iMin,iMax
43	adcroft	1.8	c(i,j)=-deltaTXrecip_hFac(i,j,1,bi,bj)recip_drF(1)
44			& KappaRX(i,j,2)recip_drC(2)
45	adcroft	1.1	b(i,j)=1.-c(i,j)
46			bet(i,j)=0.
47			IF (b(i,j).NE.0.) bet(i,j)=1. / b(i,j)
48	adcroft	1.8	gXnm1(i,j,1,bi,bj)=gXnm1(i,j,1,bi,bj)*bet(i,j)
49	adcroft	1.1	ENDDO
50			ENDDO
51			ENDIF
52	cnh	1.5	C-- Middle of forward sweep
53	cnh	1.6	IF (Nr.GT.2) THEN
54			DO k=2,Nr-1
55	adcroft	1.1	DO j=jMin,jMax
56			DO i=iMin,iMax
57			ckm1(i,j)=c(i,j)
58	adcroft	1.8	a(i,j)=-deltaTXrecip_hFac(i,j,k,bi,bj)recip_drF(k)
59			& KappaRX(i,j, k )recip_drC( k )
60			c(i,j)=-deltaTXrecip_hFac(i,j,k,bi,bj)recip_drF(k)
61			& KappaRX(i,j,k+1)recip_drC(k+1)
62	adcroft	1.1	b(i,j)=1.-c(i,j)-a(i,j)
63			ENDDO
64			ENDDO
65			DO j=jMin,jMax
66			DO i=iMin,iMax
67			gam(i,j,k)=ckm1(i,j)*bet(i,j)
68			bet(i,j)=b(i,j)-a(i,j)*gam(i,j,k)
69			IF (bet(i,j).NE.0.) bet(i,j)=1. / bet(i,j)
70	adcroft	1.8	gXnm1(i,j,k,bi,bj)=(gXnm1(i,j,k,bi,bj)
71			& -a(i,j)gXnm1(i,j,k-1,bi,bj))bet(i,j)
72	adcroft	1.1	ENDDO
73			ENDDO
74			ENDDO
75			ENDIF
76	cnh	1.6	IF (Nr.GT.1) THEN
77	cnh	1.5	C-- End of forward sweep (bottom level)
78	adcroft	1.1	DO j=jMin,jMax
79			DO i=iMin,iMax
80			ckm1(i,j)=c(i,j)
81	adcroft	1.8	a(i,j)=-deltaTXrecip_hFac(i,j,Nr,bi,bj)recip_drF(Nr)
82			& KappaRX(i,j, Nr )recip_drC( Nr )
83	adcroft	1.1	b(i,j)=1.-a(i,j)
84			ENDDO
85			ENDDO
86			DO j=jMin,jMax
87			DO i=iMin,iMax
88	cnh	1.6	gam(i,j,Nr)=ckm1(i,j)*bet(i,j)
89			bet(i,j)=b(i,j)-a(i,j)*gam(i,j,Nr)
90	adcroft	1.1	IF (bet(i,j).NE.0.) bet(i,j)=1. / bet(i,j)
91	adcroft	1.8	gXnm1(i,j,Nr,bi,bj)=(gXnm1(i,j,Nr,bi,bj)
92			& -a(i,j)gXnm1(i,j,Nr-1,bi,bj))bet(i,j)
93	adcroft	1.1	ENDDO
94			ENDDO
95	cnh	1.5	C-- Backward sweep
96	cnh	1.6	DO k=Nr-1,1,-1
97	adcroft	1.1	DO j=jMin,jMax
98			DO i=iMin,iMax
99	adcroft	1.8	gXnm1(i,j,k,bi,bj)=gXnm1(i,j,k,bi,bj)
100			& -gam(i,j,k+1)*gXnm1(i,j,k+1,bi,bj)
101	adcroft	1.1	ENDDO
102			ENDDO
103			ENDDO
104			ENDIF
105
106			RETURN
107			END