Home

Context Navigation

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

Last change on this file since 4180 was 4180, checked in by scharf, 5 years ago
removed comments in 'Former revisions' section that are older than 01.01.2019
Property svn:keywords set to `Id`
File size: 7.6 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-2019 Leibniz Universitaet Hannover
18	!------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! -----------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: sor.f90 4180 2019-08-21 14:37:54Z scharf $
27	! Rename variables in mesoscale-offline nesting mode
28	!
29	!
30	! Description:
31	! ------------
32	!> Solve the Poisson-equation with the SOR-Red/Black-scheme.
33	!------------------------------------------------------------------------------!
34	SUBROUTINE sor( d, ddzu, ddzw, p )
35
36	USE arrays_3d, &
37	ONLY: rho_air, rho_air_zw
38
39	USE grid_variables, &
40	ONLY: ddx2, ddy2
41
42	USE indices, &
43	ONLY: nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nz, nzb, nzt
44
45	USE kinds
46
47	USE control_parameters, &
48	ONLY: bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r, &
49	bc_dirichlet_s, bc_lr_cyc, bc_ns_cyc, bc_radiation_l, &
50	bc_radiation_n, bc_radiation_r, bc_radiation_s, ibc_p_b, &
51	ibc_p_t, n_sor, omega_sor
52
53	IMPLICIT NONE
54
55	INTEGER(iwp) :: i !<
56	INTEGER(iwp) :: j !<
57	INTEGER(iwp) :: k !<
58	INTEGER(iwp) :: n !<
59	INTEGER(iwp) :: nxl1 !<
60	INTEGER(iwp) :: nxl2 !<
61	INTEGER(iwp) :: nys1 !<
62	INTEGER(iwp) :: nys2 !<
63
64	REAL(wp) :: ddzu(1:nz+1) !<
65	REAL(wp) :: ddzw(1:nzt+1) !<
66
67	REAL(wp) :: d(nzb+1:nzt,nys:nyn,nxl:nxr) !<
68	REAL(wp) :: p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) !<
69
70	REAL(wp), DIMENSION(:), ALLOCATABLE :: f1 !<
71	REAL(wp), DIMENSION(:), ALLOCATABLE :: f2 !<
72	REAL(wp), DIMENSION(:), ALLOCATABLE :: f3 !<
73
74	ALLOCATE( f1(1:nz), f2(1:nz), f3(1:nz) )
75
76	!
77	!-- Compute pre-factors.
78	DO k = 1, nz
79	f2(k) = ddzu(k+1) * ddzw(k) * rho_air_zw(k)
80	f3(k) = ddzu(k) * ddzw(k) * rho_air_zw(k-1)
81	f1(k) = 2.0_wp * ( ddx2 + ddy2 ) * rho_air(k) + f2(k) + f3(k)
82	ENDDO
83
84	!
85	!-- Limits for RED- and BLACK-part.
86	IF ( MOD( nxl , 2 ) == 0 ) THEN
87	nxl1 = nxl
88	nxl2 = nxl + 1
89	ELSE
90	nxl1 = nxl + 1
91	nxl2 = nxl
92	ENDIF
93	IF ( MOD( nys , 2 ) == 0 ) THEN
94	nys1 = nys
95	nys2 = nys + 1
96	ELSE
97	nys1 = nys + 1
98	nys2 = nys
99	ENDIF
100
101	DO n = 1, n_sor
102
103	!
104	!-- RED-part
105	DO i = nxl1, nxr, 2
106	DO j = nys2, nyn, 2
107	DO k = nzb+1, nzt
108	p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * ( &
109	rho_air(k) * ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) + &
110	rho_air(k) * ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) + &
111	f2(k) * p(k+1,j,i) + &
112	f3(k) * p(k-1,j,i) - &
113	d(k,j,i) - &
114	f1(k) * p(k,j,i) )
115	ENDDO
116	ENDDO
117	ENDDO
118
119	DO i = nxl2, 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	rho_air(k) * ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) + &
124	rho_air(k) * 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	!
134	!-- Exchange of boundary values for p.
135	CALL exchange_horiz( p, nbgp )
136
137	!
138	!-- Horizontal (Neumann) boundary conditions in case of non-cyclic boundaries
139	IF ( .NOT. bc_lr_cyc ) THEN
140	IF ( bc_dirichlet_l .OR. bc_radiation_l ) p(:,:,nxl-1) = p(:,:,nxl)
141	IF ( bc_dirichlet_r .OR. bc_radiation_r ) p(:,:,nxr+1) = p(:,:,nxr)
142	ENDIF
143	IF ( .NOT. bc_ns_cyc ) THEN
144	IF ( bc_dirichlet_n .OR. bc_radiation_n ) p(:,nyn+1,:) = p(:,nyn,:)
145	IF ( bc_dirichlet_s .OR. bc_radiation_s ) p(:,nys-1,:) = p(:,nys,:)
146	ENDIF
147
148	!
149	!-- BLACK-part
150	DO i = nxl1, nxr, 2
151	DO j = nys1, nyn, 2
152	DO k = nzb+1, nzt
153	p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * ( &
154	rho_air(k) * ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) + &
155	rho_air(k) * ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) + &
156	f2(k) * p(k+1,j,i) + &
157	f3(k) * p(k-1,j,i) - &
158	d(k,j,i) - &
159	f1(k) * p(k,j,i) )
160	ENDDO
161	ENDDO
162	ENDDO
163
164	DO i = nxl2, nxr, 2
165	DO j = nys2, nyn, 2
166	DO k = nzb+1, nzt
167	p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * ( &
168	rho_air(k) * ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) + &
169	rho_air(k) * ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) + &
170	f2(k) * p(k+1,j,i) + &
171	f3(k) * p(k-1,j,i) - &
172	d(k,j,i) - &
173	f1(k) * p(k,j,i) )
174	ENDDO
175	ENDDO
176	ENDDO
177
178	!
179	!-- Exchange of boundary values for p.
180	CALL exchange_horiz( p, nbgp )
181
182	!
183	!-- Boundary conditions top/bottom.
184	!-- Bottom boundary
185	IF ( ibc_p_b == 1 ) THEN ! Neumann
186	p(nzb,:,:) = p(nzb+1,:,:)
187	ELSE ! Dirichlet
188	p(nzb,:,:) = 0.0_wp
189	ENDIF
190
191	!
192	!-- Top boundary
193	IF ( ibc_p_t == 1 ) THEN ! Neumann
194	p(nzt+1,:,:) = p(nzt,:,:)
195	ELSE ! Dirichlet
196	p(nzt+1,:,:) = 0.0_wp
197	ENDIF
198
199	!
200	!-- Horizontal (Neumann) boundary conditions in case of non-cyclic boundaries
201	IF ( .NOT. bc_lr_cyc ) THEN
202	IF ( bc_dirichlet_l .OR. bc_radiation_l ) p(:,:,nxl-1) = p(:,:,nxl)
203	IF ( bc_dirichlet_r .OR. bc_radiation_r ) p(:,:,nxr+1) = p(:,:,nxr)
204	ENDIF
205	IF ( .NOT. bc_ns_cyc ) THEN
206	IF ( bc_dirichlet_n .OR. bc_radiation_n ) p(:,nyn+1,:) = p(:,nyn,:)
207	IF ( bc_dirichlet_s .OR. bc_radiation_s ) p(:,nys-1,:) = p(:,nys,:)
208	ENDIF
209
210
211	ENDDO
212
213	DEALLOCATE( f1, f2, f3 )
214
215	END SUBROUTINE sor

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |