Home

Context Navigation

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

Last change on this file since 1947 was 1818, checked in by maronga, 9 years ago
last commit documented / copyright update
Property svn:keywords set to `Id`
File size: 8.0 KB

Line
1	!> @file sor.f90
2	!--------------------------------------------------------------------------------!
3	! This file is part of PALM.
4	!
5	! PALM is free software: you can redistribute it and/or modify it under the terms
6	! of the GNU General Public License as published by the Free Software Foundation,
7	! either version 3 of the License, or (at your option) any later version.
8	!
9	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
10	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
11	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
12	!
13	! You should have received a copy of the GNU General Public License along with
14	! PALM. If not, see <http://www.gnu.org/licenses/>.
15	!
16	! Copyright 1997-2016 Leibniz Universitaet Hannover
17	!--------------------------------------------------------------------------------!
18	!
19	! Current revisions:
20	! -----------------
21	!
22	!
23	! Former revisions:
24	! -----------------
25	! $Id: sor.f90 1818 2016-04-06 15:53:27Z suehring $
26	!
27	! 1762 2016-02-25 12:31:13Z hellstea
28	! Introduction of nested domain feature
29	!
30	! 1682 2015-10-07 23:56:08Z knoop
31	! Code annotations made doxygen readable
32	!
33	! 1353 2014-04-08 15:21:23Z heinze
34	! REAL constants provided with KIND-attribute
35	!
36	! 1320 2014-03-20 08:40:49Z raasch
37	! ONLY-attribute added to USE-statements,
38	! kind-parameters added to all INTEGER and REAL declaration statements,
39	! kinds are defined in new module kinds,
40	! old module precision_kind is removed,
41	! revision history before 2012 removed,
42	! comment fields (!:) to be used for variable explanations added to
43	! all variable declaration statements
44	!
45	! 1036 2012-10-22 13:43:42Z raasch
46	! code put under GPL (PALM 3.9)
47	!
48	! Revision 1.1 1997/08/11 06:25:56 raasch
49	! Initial revision
50	!
51	!
52	! Description:
53	! ------------
54	!> Solve the Poisson-equation with the SOR-Red/Black-scheme.
55	!------------------------------------------------------------------------------!
56	SUBROUTINE sor( d, ddzu, ddzw, p )
57
58
59	USE grid_variables, &
60	ONLY: ddx2, ddy2
61
62	USE indices, &
63	ONLY: nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nz, nzb, nzt
64
65	USE kinds
66
67	USE control_parameters, &
68	ONLY: bc_lr_cyc, bc_ns_cyc, ibc_p_b, ibc_p_t, inflow_l, inflow_n, &
69	inflow_r, inflow_s, nest_bound_l, nest_bound_n, nest_bound_r, &
70	nest_bound_s, n_sor, omega_sor, outflow_l, outflow_n, &
71	outflow_r, outflow_s
72
73	IMPLICIT NONE
74
75	INTEGER(iwp) :: i !<
76	INTEGER(iwp) :: j !<
77	INTEGER(iwp) :: k !<
78	INTEGER(iwp) :: n !<
79	INTEGER(iwp) :: nxl1 !<
80	INTEGER(iwp) :: nxl2 !<
81	INTEGER(iwp) :: nys1 !<
82	INTEGER(iwp) :: nys2 !<
83
84	REAL(wp) :: ddzu(1:nz+1) !<
85	REAL(wp) :: ddzw(1:nzt+1) !<
86
87	REAL(wp) :: d(nzb+1:nzt,nys:nyn,nxl:nxr) !<
88	REAL(wp) :: p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) !<
89
90	REAL(wp), DIMENSION(:), ALLOCATABLE :: f1 !<
91	REAL(wp), DIMENSION(:), ALLOCATABLE :: f2 !<
92	REAL(wp), DIMENSION(:), ALLOCATABLE :: f3 !<
93
94	ALLOCATE( f1(1:nz), f2(1:nz), f3(1:nz) )
95
96	!
97	!-- Compute pre-factors.
98	DO k = 1, nz
99	f2(k) = ddzu(k+1) * ddzw(k)
100	f3(k) = ddzu(k) * ddzw(k)
101	f1(k) = 2.0_wp * ( ddx2 + ddy2 ) + f2(k) + f3(k)
102	ENDDO
103
104	!
105	!-- Limits for RED- and BLACK-part.
106	IF ( MOD( nxl , 2 ) == 0 ) THEN
107	nxl1 = nxl
108	nxl2 = nxl + 1
109	ELSE
110	nxl1 = nxl + 1
111	nxl2 = nxl
112	ENDIF
113	IF ( MOD( nys , 2 ) == 0 ) THEN
114	nys1 = nys
115	nys2 = nys + 1
116	ELSE
117	nys1 = nys + 1
118	nys2 = nys
119	ENDIF
120
121	DO n = 1, n_sor
122
123	!
124	!-- RED-part
125	DO i = nxl1, nxr, 2
126	DO j = nys2, nyn, 2
127	DO k = nzb+1, nzt
128	p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * ( &
129	ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) + &
130	ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) + &
131	f2(k) * p(k+1,j,i) + &
132	f3(k) * p(k-1,j,i) - &
133	d(k,j,i) - &
134	f1(k) * p(k,j,i) )
135	ENDDO
136	ENDDO
137	ENDDO
138
139	DO i = nxl2, nxr, 2
140	DO j = nys1, nyn, 2
141	DO k = nzb+1, nzt
142	p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * ( &
143	ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) + &
144	ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) + &
145	f2(k) * p(k+1,j,i) + &
146	f3(k) * p(k-1,j,i) - &
147	d(k,j,i) - &
148	f1(k) * p(k,j,i) )
149	ENDDO
150	ENDDO
151	ENDDO
152
153	!
154	!-- Exchange of boundary values for p.
155	CALL exchange_horiz( p, nbgp )
156
157	!
158	!-- Horizontal (Neumann) boundary conditions in case of non-cyclic boundaries
159	IF ( .NOT. bc_lr_cyc ) THEN
160	IF ( inflow_l .OR. outflow_l .OR. nest_bound_l ) p(:,:,nxl-1) = p(:,:,nxl)
161	IF ( inflow_r .OR. outflow_r .OR. nest_bound_r ) p(:,:,nxr+1) = p(:,:,nxr)
162	ENDIF
163	IF ( .NOT. bc_ns_cyc ) THEN
164	IF ( inflow_n .OR. outflow_n .OR. nest_bound_n ) p(:,nyn+1,:) = p(:,nyn,:)
165	IF ( inflow_s .OR. outflow_s .OR. nest_bound_s ) p(:,nys-1,:) = p(:,nys,:)
166	ENDIF
167
168	!
169	!-- BLACK-part
170	DO i = nxl1, nxr, 2
171	DO j = nys1, nyn, 2
172	DO k = nzb+1, nzt
173	p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * ( &
174	ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) + &
175	ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) + &
176	f2(k) * p(k+1,j,i) + &
177	f3(k) * p(k-1,j,i) - &
178	d(k,j,i) - &
179	f1(k) * p(k,j,i) )
180	ENDDO
181	ENDDO
182	ENDDO
183
184	DO i = nxl2, nxr, 2
185	DO j = nys2, nyn, 2
186	DO k = nzb+1, nzt
187	p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * ( &
188	ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) + &
189	ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) + &
190	f2(k) * p(k+1,j,i) + &
191	f3(k) * p(k-1,j,i) - &
192	d(k,j,i) - &
193	f1(k) * p(k,j,i) )
194	ENDDO
195	ENDDO
196	ENDDO
197
198	!
199	!-- Exchange of boundary values for p.
200	CALL exchange_horiz( p, nbgp )
201
202	!
203	!-- Boundary conditions top/bottom.
204	!-- Bottom boundary
205	IF ( ibc_p_b == 1 ) THEN ! Neumann
206	p(nzb,:,:) = p(nzb+1,:,:)
207	ELSE ! Dirichlet
208	p(nzb,:,:) = 0.0_wp
209	ENDIF
210
211	!
212	!-- Top boundary
213	IF ( ibc_p_t == 1 ) THEN ! Neumann
214	p(nzt+1,:,:) = p(nzt,:,:)
215	ELSE ! Dirichlet
216	p(nzt+1,:,:) = 0.0_wp
217	ENDIF
218
219	!
220	!-- Horizontal (Neumann) boundary conditions in case of non-cyclic boundaries
221	IF ( .NOT. bc_lr_cyc ) THEN
222	IF ( inflow_l .OR. outflow_l .OR. nest_bound_l ) p(:,:,nxl-1) = p(:,:,nxl)
223	IF ( inflow_r .OR. outflow_r .OR. nest_bound_r ) p(:,:,nxr+1) = p(:,:,nxr)
224	ENDIF
225	IF ( .NOT. bc_ns_cyc ) THEN
226	IF ( inflow_n .OR. outflow_n .OR. nest_bound_n ) p(:,nyn+1,:) = p(:,nyn,:)
227	IF ( inflow_s .OR. outflow_s .OR. nest_bound_s ) p(:,nys-1,:) = p(:,nys,:)
228	ENDIF
229
230
231	ENDDO
232
233	DEALLOCATE( f1, f2, f3 )
234
235	END SUBROUTINE sor

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |