/[MITgcm]/MITgcm_contrib/gael/pkg/smooth2/smooth_impldiff.F
ViewVC logotype

Annotation of /MITgcm_contrib/gael/pkg/smooth2/smooth_impldiff.F

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

Revision 1.1 - (hide annotations) (download)
Sun Oct 25 21:24:04 2009 UTC (15 years, 8 months ago) by gforget
Branch: MAIN
CVS Tags: HEAD 
Renovated pkg/smooth. Ready for MITgcm/pkg check in?
1 gforget 1.1 C $Header: /u/gcmpack/MITgcm_contrib/gael/pkg/smooth2/smooth_impldiff.F,v 1.1 2009/10/24 23:27:24 gforget Exp $
2     C $Name: $
3    
4     #include "PACKAGES_CONFIG.h"
5     #include "CPP_OPTIONS.h"
6    
7     SUBROUTINE SMOOTH_IMPLDIFF( bi, bj, iMin, iMax, jMin, jMax,
8     I deltaTX, KappaRX, recip_hFac,
9     U gXNm1,
10     I myThid )
11    
12     C *==========================================================*
13     C | SUBROUTINE smooth_impldiff
14     C | o Copy of model/src/impldiff
15     C | (simplified and with specified time step)
16     C *==========================================================*
17    
18    
19     C !USES:
20     IMPLICIT NONE
21     C == Global data ==
22     #include "SIZE.h"
23     #include "DYNVARS.h"
24     #include "EEPARAMS.h"
25     #include "PARAMS.h"
26     #include "GRID.h"
27     #ifdef ALLOW_GENERIC_ADVDIFF
28     #include "GAD.h"
29     #endif
30     #ifdef ALLOW_LONGSTEP
31     #include "LONGSTEP_PARAMS.h"
32     #endif
33     #ifdef ALLOW_PTRACERS
34     #include "PTRACERS_SIZE.h"
35     #include "PTRACERS_PARAMS.h"
36     #endif
37     #ifdef ALLOW_AUTODIFF_TAMC
38     #include "tamc_keys.h"
39     #endif
40    
41     C !INPUT/OUTPUT PARAMETERS:
42     C == Routine Arguments ==
43     INTEGER bi,bj,iMin,iMax,jMin,jMax
44     _RL KappaRX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
45     _RS recip_hFac(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
46     _RL gXnm1(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
47     INTEGER myThid
48    
49     C !LOCAL VARIABLES:
50     C == Local variables ==
51     INTEGER i,j,k
52     _RL deltaTX
53     _RL gYnm1(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
54     _RL a(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
55     _RL b(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
56     _RL c(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
57     _RL bet(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
58     _RL gam(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
59     CEOP
60    
61     c deltaTX = 1. _d 0
62    
63     C-- Initialise
64     DO k=1,Nr
65     DO j=jMin,jMax
66     DO i=iMin,iMax
67     gYNm1(i,j,k,bi,bj) = 0. _d 0
68     ENDDO
69     ENDDO
70     ENDDO
71    
72     C-- Old aLower
73     DO j=jMin,jMax
74     DO i=iMin,iMax
75     a(i,j,1) = 0. _d 0
76     ENDDO
77     ENDDO
78     DO k=2,Nr
79     DO j=jMin,jMax
80     DO i=iMin,iMax
81     a(i,j,k) = -deltaTX*recip_hFac(i,j,k,bi,bj)*recip_drF(k)
82     & *recip_deepFac2C(k)*recip_rhoFacC(k)
83     & *KappaRX(i,j, k )*recip_drC( k )
84     & *deepFac2F(k)*rhoFacF(k)
85     IF (recip_hFac(i,j,k-1,bi,bj).EQ.0.) a(i,j,k)=0.
86     ENDDO
87     ENDDO
88     ENDDO
89    
90     C-- Old aUpper
91     DO k=1,Nr-1
92     DO j=jMin,jMax
93     DO i=iMin,iMax
94     c(i,j,k) = -deltaTX*recip_hFac(i,j,k,bi,bj)*recip_drF(k)
95     & *recip_deepFac2C(k)*recip_rhoFacC(k)
96     & *KappaRX(i,j,k+1)*recip_drC(k+1)
97     & *deepFac2F(k+1)*rhoFacF(k+1)
98     IF (recip_hFac(i,j,k+1,bi,bj).EQ.0.) c(i,j,k)=0.
99     ENDDO
100     ENDDO
101     ENDDO
102     DO j=jMin,jMax
103     DO i=iMin,iMax
104     c(i,j,Nr) = 0. _d 0
105     ENDDO
106     ENDDO
107    
108     C-- Old aCenter
109     DO k=1,Nr
110     DO j=jMin,jMax
111     DO i=iMin,iMax
112     b(i,j,k) = 1. _d 0 - c(i,j,k) - a(i,j,k)
113     ENDDO
114     ENDDO
115     ENDDO
116    
117     C-- Old and new gam, bet are the same
118     DO k=1,Nr
119     DO j=jMin,jMax
120     DO i=iMin,iMax
121     bet(i,j,k) = 1. _d 0
122     gam(i,j,k) = 0. _d 0
123     ENDDO
124     ENDDO
125     ENDDO
126    
127     C-- Only need do anything if Nr>1
128     IF (Nr.GT.1) THEN
129    
130     k = 1
131     C-- Beginning of forward sweep (top level)
132     DO j=jMin,jMax
133     DO i=iMin,iMax
134     IF (b(i,j,1).NE.0.) bet(i,j,1) = 1. _d 0 / b(i,j,1)
135     ENDDO
136     ENDDO
137    
138     ENDIF
139    
140     C-- Middle of forward sweep
141     IF (Nr.GE.2) THEN
142    
143     CADJ loop = sequential
144     DO k=2,Nr
145    
146     DO j=jMin,jMax
147     DO i=iMin,iMax
148     gam(i,j,k) = c(i,j,k-1)*bet(i,j,k-1)
149     IF ( ( b(i,j,k) - a(i,j,k)*gam(i,j,k) ) .NE. 0.)
150     & bet(i,j,k) = 1. _d 0 / ( b(i,j,k) - a(i,j,k)*gam(i,j,k) )
151     ENDDO
152     ENDDO
153    
154     ENDDO
155    
156     ENDIF
157    
158    
159     DO j=jMin,jMax
160     DO i=iMin,iMax
161     gYNm1(i,j,1,bi,bj) = gXNm1(i,j,1,bi,bj)*bet(i,j,1)
162     ENDDO
163     ENDDO
164     DO k=2,Nr
165     DO j=jMin,jMax
166     DO i=iMin,iMax
167     gYnm1(i,j,k,bi,bj) = bet(i,j,k)*
168     & (gXnm1(i,j,k,bi,bj) - a(i,j,k)*gYnm1(i,j,k-1,bi,bj))
169     ENDDO
170     ENDDO
171     ENDDO
172    
173    
174     C-- Backward sweep
175     CADJ loop = sequential
176     DO k=Nr-1,1,-1
177     DO j=jMin,jMax
178     DO i=iMin,iMax
179     gYnm1(i,j,k,bi,bj)=gYnm1(i,j,k,bi,bj)
180     & -gam(i,j,k+1)*gYnm1(i,j,k+1,bi,bj)
181     ENDDO
182     ENDDO
183     ENDDO
184    
185     DO k=1,Nr
186     DO j=jMin,jMax
187     DO i=iMin,iMax
188     gXnm1(i,j,k,bi,bj)=gYnm1(i,j,k,bi,bj)
189     ENDDO
190     ENDDO
191     ENDDO
192    
193     RETURN
194     END
  ViewVC Help
Powered by ViewVC 1.1.22