Home

Context Navigation

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

Last change on this file since 4516 was 4457, checked in by raasch, 5 years ago
ghost point exchange modularized, bugfix for wrong 2d-exchange
Property svn:keywords set to `Id`
File size: 7.9 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |