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

Last change on this file since 4651 was 4591, checked in by raasch, 4 years ago

files re-formatted to follow the PALM coding standard

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