model/src/impldiff.F

C $Header: /u/gcmpack/models/MITgcmUV/model/src/impldiff.F,v 1.15 2001/02/04 14:38:47 cnh Exp $
C $Name: checkpoint38 $

#include "CPP_OPTIONS.h"

C     |==========================================================|
C     | S/R IMPLDIFF                                             |
C     | o Solve implicit diffusion equation for vertical         |
C     |   diffusivity.                                           |
C     | o Recoded from 2d intermediate fields to 3d to reduce    |
C     |   TAMC storage                                           |
C     | o Fixed missing masks for fields a(), c()                |
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)
      _RL gYnm1(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,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)

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*recip_hFac(i,j,k,bi,bj)*recip_drF(k)
     &               *KappaRX(i,j, k )*recip_drC( k )
        ENDDO
       ENDDO
      ENDDO

C--   Old aUpper
      DO k=1,Nr-1
       DO j=jMin,jMax
        DO i=iMin,iMax
          c(i,j,k) = -deltaTX*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) = 0. _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.GT.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

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