Home

Context Navigation

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

Last change on this file since 4397 was 4360, checked in by suehring, 5 years ago
Bugfix in output of time-averaged plant-canopy quanities; Output of plant-canopy data only where tall canopy is defined; land-surface model: fix wrong location strings; tests: update urban test case; all source code files: copyright update
Property svn:keywords set to `Id`
File size: 7.7 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 4360 2020-01-07 11:25:50Z scharf $
27	! Corrected "Former revisions" section
28	!
29	! 3655 2019-01-07 16:51:22Z knoop
30	! Rename variables in mesoscale-offline nesting mode
31	!
32	! Revision 1.1 1997/08/11 06:25:56 raasch
33	! Initial revision
34	!
35	!
36	! Description:
37	! ------------
38	!> Solve the Poisson-equation with the SOR-Red/Black-scheme.
39	!------------------------------------------------------------------------------!
40	SUBROUTINE sor( d, ddzu, ddzw, p )
41
42	USE arrays_3d, &
43	ONLY: rho_air, rho_air_zw
44
45	USE grid_variables, &
46	ONLY: ddx2, ddy2
47
48	USE indices, &
49	ONLY: nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nz, nzb, nzt
50
51	USE kinds
52
53	USE control_parameters, &
54	ONLY: bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r, &
55	bc_dirichlet_s, bc_lr_cyc, bc_ns_cyc, bc_radiation_l, &
56	bc_radiation_n, bc_radiation_r, bc_radiation_s, ibc_p_b, &
57	ibc_p_t, n_sor, omega_sor
58
59	IMPLICIT NONE
60
61	INTEGER(iwp) :: i !<
62	INTEGER(iwp) :: j !<
63	INTEGER(iwp) :: k !<
64	INTEGER(iwp) :: n !<
65	INTEGER(iwp) :: nxl1 !<
66	INTEGER(iwp) :: nxl2 !<
67	INTEGER(iwp) :: nys1 !<
68	INTEGER(iwp) :: nys2 !<
69
70	REAL(wp) :: ddzu(1:nz+1) !<
71	REAL(wp) :: ddzw(1:nzt+1) !<
72
73	REAL(wp) :: d(nzb+1:nzt,nys:nyn,nxl:nxr) !<
74	REAL(wp) :: p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) !<
75
76	REAL(wp), DIMENSION(:), ALLOCATABLE :: f1 !<
77	REAL(wp), DIMENSION(:), ALLOCATABLE :: f2 !<
78	REAL(wp), DIMENSION(:), ALLOCATABLE :: f3 !<
79
80	ALLOCATE( f1(1:nz), f2(1:nz), f3(1:nz) )
81
82	!
83	!-- Compute pre-factors.
84	DO k = 1, nz
85	f2(k) = ddzu(k+1) * ddzw(k) * rho_air_zw(k)
86	f3(k) = ddzu(k) * ddzw(k) * rho_air_zw(k-1)
87	f1(k) = 2.0_wp * ( ddx2 + ddy2 ) * rho_air(k) + f2(k) + f3(k)
88	ENDDO
89
90	!
91	!-- Limits for RED- and BLACK-part.
92	IF ( MOD( nxl , 2 ) == 0 ) THEN
93	nxl1 = nxl
94	nxl2 = nxl + 1
95	ELSE
96	nxl1 = nxl + 1
97	nxl2 = nxl
98	ENDIF
99	IF ( MOD( nys , 2 ) == 0 ) THEN
100	nys1 = nys
101	nys2 = nys + 1
102	ELSE
103	nys1 = nys + 1
104	nys2 = nys
105	ENDIF
106
107	DO n = 1, n_sor
108
109	!
110	!-- RED-part
111	DO i = nxl1, nxr, 2
112	DO j = nys2, nyn, 2
113	DO k = nzb+1, nzt
114	p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * ( &
115	rho_air(k) * ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) + &
116	rho_air(k) * ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) + &
117	f2(k) * p(k+1,j,i) + &
118	f3(k) * p(k-1,j,i) - &
119	d(k,j,i) - &
120	f1(k) * p(k,j,i) )
121	ENDDO
122	ENDDO
123	ENDDO
124
125	DO i = nxl2, nxr, 2
126	DO j = nys1, nyn, 2
127	DO k = nzb+1, nzt
128	p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * ( &
129	rho_air(k) * ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) + &
130	rho_air(k) * 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	!
140	!-- Exchange of boundary values for p.
141	CALL exchange_horiz( p, nbgp )
142
143	!
144	!-- Horizontal (Neumann) boundary conditions in case of non-cyclic boundaries
145	IF ( .NOT. bc_lr_cyc ) THEN
146	IF ( bc_dirichlet_l .OR. bc_radiation_l ) p(:,:,nxl-1) = p(:,:,nxl)
147	IF ( bc_dirichlet_r .OR. bc_radiation_r ) p(:,:,nxr+1) = p(:,:,nxr)
148	ENDIF
149	IF ( .NOT. bc_ns_cyc ) THEN
150	IF ( bc_dirichlet_n .OR. bc_radiation_n ) p(:,nyn+1,:) = p(:,nyn,:)
151	IF ( bc_dirichlet_s .OR. bc_radiation_s ) p(:,nys-1,:) = p(:,nys,:)
152	ENDIF
153
154	!
155	!-- BLACK-part
156	DO i = nxl1, nxr, 2
157	DO j = nys1, nyn, 2
158	DO k = nzb+1, nzt
159	p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * ( &
160	rho_air(k) * ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) + &
161	rho_air(k) * ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) + &
162	f2(k) * p(k+1,j,i) + &
163	f3(k) * p(k-1,j,i) - &
164	d(k,j,i) - &
165	f1(k) * p(k,j,i) )
166	ENDDO
167	ENDDO
168	ENDDO
169
170	DO i = nxl2, nxr, 2
171	DO j = nys2, nyn, 2
172	DO k = nzb+1, nzt
173	p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * ( &
174	rho_air(k) * ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) + &
175	rho_air(k) * 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	!
185	!-- Exchange of boundary values for p.
186	CALL exchange_horiz( p, nbgp )
187
188	!
189	!-- Boundary conditions top/bottom.
190	!-- Bottom boundary
191	IF ( ibc_p_b == 1 ) THEN ! Neumann
192	p(nzb,:,:) = p(nzb+1,:,:)
193	ELSE ! Dirichlet
194	p(nzb,:,:) = 0.0_wp
195	ENDIF
196
197	!
198	!-- Top boundary
199	IF ( ibc_p_t == 1 ) THEN ! Neumann
200	p(nzt+1,:,:) = p(nzt,:,:)
201	ELSE ! Dirichlet
202	p(nzt+1,:,:) = 0.0_wp
203	ENDIF
204
205	!
206	!-- Horizontal (Neumann) boundary conditions in case of non-cyclic boundaries
207	IF ( .NOT. bc_lr_cyc ) THEN
208	IF ( bc_dirichlet_l .OR. bc_radiation_l ) p(:,:,nxl-1) = p(:,:,nxl)
209	IF ( bc_dirichlet_r .OR. bc_radiation_r ) p(:,:,nxr+1) = p(:,:,nxr)
210	ENDIF
211	IF ( .NOT. bc_ns_cyc ) THEN
212	IF ( bc_dirichlet_n .OR. bc_radiation_n ) p(:,nyn+1,:) = p(:,nyn,:)
213	IF ( bc_dirichlet_s .OR. bc_radiation_s ) p(:,nys-1,:) = p(:,nys,:)
214	ENDIF
215
216
217	ENDDO
218
219	DEALLOCATE( f1, f2, f3 )
220
221	END SUBROUTINE sor

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |