Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

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

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

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