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

Last change on this file since 698 was 668, checked in by suehring, 13 years ago

last commit documented

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