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

Last change on this file since 504 was 484, checked in by raasch, 15 years ago

typo in file headers removed

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