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

Last change on this file since 707 was 707, checked in by raasch, 10 years ago

New:
---

In case of multigrid method, on coarse grid levels, gathered data are
identically processed on all PEs (before, on PE0 only), so that the subsequent
scattering of data is not neccessary any more. (modules, init_pegrid, poismg)

Changed:


Calculation of weighted average of p is now handled in the same way
regardless of the number of ghost layers (advection scheme). (pres)

multigrid and sor method are using p_loc for iterative
advancements of pressure. p_sub removed. (init_3d_model, modules, poismg, pres, sor)

bc_lr and bc_ns replaced by bc_lr_dirrad, bc_lr_raddir, bc_ns_dirrad, bc_ns_raddir
for speed optimization. (calc_spectra, check_parameters, exchange_horiz,
exchange_horiz_2d, header, init_3d_model, init_grid, init_pegrid, modules,
poismg, pres, sor, time_integration, timestep)

grid_level directly used as index for MPI data type arrays. (exchange_horiz,
poismg)

initial assignments of zero to array p for iterative solvers only (init_3d_model)

Errors:


localsum calculation modified for proper OpenMP reduction. (pres)

Bugfix: bottom (nzb) and top (nzt+1) boundary conditions set in routines
resid and restrict. They were missed before, which may have led to
unpredictable results. (poismg)

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