source: palm/tags/release-3.1b/SOURCE/sor.f90 @ 710

Last change on this file since 710 was 4, checked in by raasch, 17 years ago

Id keyword set as property for all *.f90 files

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