source: palm/trunk/SOURCE/sor.f90 @ 1

Last change on this file since 1 was 1, checked in by raasch, 14 years ago

Initial repository layout and content

File size: 6.2 KB
Line 
1 SUBROUTINE sor( d, ddzu, ddzw, p )
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Log: sor.f90,v $
11! Revision 1.9  2005/03/26 21:02:23  raasch
12! Implementation of non-cyclic (Neumann) horizontal boundary conditions,
13! dx2,dy2 replaced by ddx2,ddy2
14!
15! Revision 1.8  2001/03/30 07:52:44  raasch
16! Arrays f1-f3 changed from dummy arguments to allocatable arrays,
17! Translation of remaining German identifiers (variables, subroutines, etc.)
18!
19! Revision 1.7  2001/01/30 21:41:25  letzel
20! All comments translated into English.
21!
22! Revision 1.6  2001/01/22 08:06:14  raasch
23! Module test_variables removed
24!
25! Revision 1.5  1998/07/06 12:33:13  raasch
26! + USE test_variables
27!
28! Revision 1.4  1997/09/12 06:30:12  raasch
29! Randbedingungen umbenannt
30!
31! Revision 1.3  1997/09/09 08:30:44  raasch
32! Kehrwerte der Gitterweiten implementiert
33!
34! Revision 1.2  1997/08/29 09:00:58  raasch
35! omega --> omega_sor
36!
37! Revision 1.1  1997/08/11 06:25:56  raasch
38! Initial revision
39!
40!
41! Description:
42! ------------
43! Solve the Poisson-equation with the SOR-Red/Black-scheme.
44!-------------------------------------------------------------------------------!
45
46    USE grid_variables
47    USE indices
48    USE pegrid
49    USE control_parameters
50
51    IMPLICIT NONE
52
53    INTEGER ::  i, j, k, n, nxl1, nxl2, nys1, nys2
54    REAL    ::  ddzu(1:nz+1), ddzw(1:nz)
55    REAL    ::  d(nzb+1:nzt,nys:nyn,nxl:nxr),         &
56                p(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
57    REAL, DIMENSION(:), ALLOCATABLE ::  f1, f2, f3
58
59    ALLOCATE( f1(1:nz), f2(1:nz), f3(1:nz) )
60
61!
62!-- Compute pre-factors.
63    DO  k = 1, nz
64         f2(k) = ddzu(k+1) * ddzw(k)
65         f3(k) = ddzu(k)   * ddzw(k)
66         f1(k) = 2.0 * ( ddx2 + ddy2 ) + f2(k) + f3(k)
67    ENDDO
68
69!
70!-- Limits for RED- and BLACK-part.
71    IF ( MOD( nxl , 2 ) == 0 )  THEN
72       nxl1 = nxl
73       nxl2 = nxl + 1
74    ELSE
75       nxl1 = nxl + 1
76       nxl2 = nxl
77    ENDIF
78    IF ( MOD( nys , 2 ) == 0 )  THEN
79       nys1 = nys
80       nys2 = nys + 1
81    ELSE
82       nys1 = nys + 1
83       nys2 = nys
84    ENDIF
85
86    DO  n = 1, n_sor
87
88!
89!--    RED-part
90       DO  i = nxl1, nxr, 2
91          DO  j = nys2, nyn, 2
92             DO  k = nzb+1, nzt
93                p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * (            &
94                               ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +    &
95                               ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +    &
96                               f2(k) * p(k+1,j,i)                 +    &
97                               f3(k) * p(k-1,j,i)                 -    &
98                               d(k,j,i)                           -    &
99                               f1(k) * p(k,j,i)           )
100             ENDDO
101          ENDDO
102       ENDDO
103
104       DO  i = nxl2, nxr, 2
105          DO  j = nys1, nyn, 2
106             DO  k = nzb+1, nzt
107                p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * (            &
108                               ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +    &
109                               ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +    &
110                               f2(k) * p(k+1,j,i)                 +    &
111                               f3(k) * p(k-1,j,i)                 -    &
112                               d(k,j,i)                           -    &
113                               f1(k) * p(k,j,i)           )
114             ENDDO
115          ENDDO
116       ENDDO
117
118!
119!--    Exchange of boundary values for p.
120       CALL exchange_horiz( p, 0, 0 )
121
122!
123!--    Horizontal (Neumann) boundary conditions in case of non-cyclic boundaries
124       IF ( bc_lr /= 'cyclic' )  THEN
125          IF ( inflow_l .OR. outflow_l )  p(:,:,nxl-1) = p(:,:,nxl)
126          IF ( inflow_r .OR. outflow_r )  p(:,:,nxr+1) = p(:,:,nxr)
127       ENDIF
128       IF ( bc_ns /= 'cyclic' )  THEN
129          IF ( inflow_n .OR. outflow_n )  p(:,nyn+1,:) = p(:,nyn,:)
130          IF ( inflow_s .OR. outflow_s )  p(:,nys-1,:) = p(:,nys,:)
131       ENDIF
132
133!
134!--    BLACK-part
135       DO  i = nxl1, nxr, 2
136          DO  j = nys1, nyn, 2
137             DO  k = nzb+1, nzt
138                p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * (            &
139                               ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +    &
140                               ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +    &
141                               f2(k) * p(k+1,j,i)                 +    &
142                               f3(k) * p(k-1,j,i)                 -    &
143                               d(k,j,i)                           -    &
144                               f1(k) * p(k,j,i)           )
145             ENDDO
146          ENDDO
147       ENDDO
148
149       DO  i = nxl2, nxr, 2
150          DO  j = nys2, nyn, 2
151             DO  k = nzb+1, nzt
152                p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * (            &
153                               ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +    &
154                               ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +    &
155                               f2(k) * p(k+1,j,i)                 +    &
156                               f3(k) * p(k-1,j,i)                 -    &
157                               d(k,j,i)                           -    &
158                               f1(k) * p(k,j,i)           )
159             ENDDO
160          ENDDO
161       ENDDO
162
163!
164!--    Exchange of boundary values for p.
165       CALL exchange_horiz( p, 0, 0 )
166
167!
168!--    Boundary conditions top/bottom.
169!--    Bottom boundary
170       IF ( ibc_p_b == 1 )  THEN
171!
172!--       Neumann
173          p(nzb,:,:) = p(nzb+1,:,:)
174       ELSE
175!
176!--       Dirichlet
177          p(nzb,:,:) = 0.0
178       ENDIF
179
180!
181!--    Top boundary
182       IF ( ibc_p_t == 1 )  THEN
183!
184!--       Neumann
185          p(nzt+1,:,:) = p(nzt,:,:)
186       ELSE
187!
188!--       Dirichlet
189          p(nzt+1,:,:) = 0.0
190       ENDIF
191
192!
193!--    Horizontal (Neumann) boundary conditions in case of non-cyclic boundaries
194       IF ( bc_lr /= 'cyclic' )  THEN
195          IF ( inflow_l .OR. outflow_l )  p(:,:,nxl-1) = p(:,:,nxl)
196          IF ( inflow_r .OR. outflow_r )  p(:,:,nxr+1) = p(:,:,nxr)
197       ENDIF
198       IF ( bc_ns /= 'cyclic' )  THEN
199          IF ( inflow_n .OR. outflow_n )  p(:,nyn+1,:) = p(:,nyn,:)
200          IF ( inflow_s .OR. outflow_s )  p(:,nys-1,:) = p(:,nys,:)
201       ENDIF
202
203    ENDDO
204
205    DEALLOCATE( f1, f2, f3 )
206
207 END SUBROUTINE sor
Note: See TracBrowser for help on using the repository browser.