/[MITgcm]/MITgcm/model/src/impldiff.F
ViewVC logotype

Contents of /MITgcm/model/src/impldiff.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.14 - (show annotations) (download)
Mon Jan 29 20:04:43 2001 UTC (23 years, 4 months ago) by heimbach
Branch: MAIN
CVS Tags: checkpoint34
Changes since 1.13: +15 -23 lines
Changed loop boundaries.

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

  ViewVC Help
Powered by ViewVC 1.1.22