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

Last change on this file since 4879 was 4828, checked in by Giersch, 4 years ago

Copyright updated to year 2021, interface pmc_sort removed to accelarate the nesting code

  • Property svn:keywords set to Id
File size: 10.9 KB
RevLine 
[1682]1!> @file sor.f90
[4591]2!--------------------------------------------------------------------------------------------------!
[2696]3! This file is part of the PALM model system.
[1036]4!
[4591]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.
[1036]8!
[4591]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.
[1036]12!
[4591]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/>.
[1036]15!
[4828]16! Copyright 1997-2021 Leibniz Universitaet Hannover
[4591]17!--------------------------------------------------------------------------------------------------!
[1036]18!
[4591]19!
[484]20! Current revisions:
[1]21! -----------------
[4591]22!
23!
[1321]24! Former revisions:
25! -----------------
26! $Id: sor.f90 4828 2021-01-05 11:21:41Z gronemeier $
[4591]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!
[4457]33! 4360 2020-01-07 11:25:50Z suehring
[4182]34! Corrected "Former revisions" section
[4591]35!
[4182]36! 3655 2019-01-07 16:51:22Z knoop
[3183]37! Rename variables in mesoscale-offline nesting mode
[1321]38!
[4182]39! Revision 1.1  1997/08/11 06:25:56  raasch
40! Initial revision
41!
[4591]42!--------------------------------------------------------------------------------------------------!
[1]43! Description:
44! ------------
[1682]45!> Solve the Poisson-equation with the SOR-Red/Black-scheme.
[4591]46!--------------------------------------------------------------------------------------------------!
[1682]47 SUBROUTINE sor( d, ddzu, ddzw, p )
[1]48
[4591]49    USE arrays_3d,                                                                                 &
50        ONLY:  rho_air,                                                                            &
51               rho_air_zw
[2037]52
[4591]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
[4457]68
[4591]69    USE exchange_horiz_mod,                                                                        &
[4457]70        ONLY:  exchange_horiz
71
[4591]72    USE grid_variables,                                                                            &
73        ONLY:  ddx2,                                                                               &
74               ddy2
[1]75
[4591]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
[1320]89
90    USE kinds
91
[1]92    IMPLICIT NONE
93
[4591]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  !<
[1]102
[4591]103    REAL(wp) ::  ddzu(1:nz+1)   !<
104    REAL(wp) ::  ddzw(1:nzt+1)  !<
[1320]105
[4591]106    REAL(wp) ::  d(nzb+1:nzt,nys:nyn,nxl:nxr)      !<
107    REAL(wp) ::  p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  !<
[1320]108
[4591]109    REAL(wp), DIMENSION(:), ALLOCATABLE ::  f1  !<
110    REAL(wp), DIMENSION(:), ALLOCATABLE ::  f2  !<
111    REAL(wp), DIMENSION(:), ALLOCATABLE ::  f3  !<
[1320]112
[1]113    ALLOCATE( f1(1:nz), f2(1:nz), f3(1:nz) )
114
115!
116!-- Compute pre-factors.
117    DO  k = 1, nz
[2037]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)
[1]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
[4591]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)               )
[1]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
[4591]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)               )
[1]168             ENDDO
169          ENDDO
170       ENDDO
171
172!
173!--    Exchange of boundary values for p.
[667]174       CALL exchange_horiz( p, nbgp )
[1]175
[4591]176
[1]177!
178!--    Horizontal (Neumann) boundary conditions in case of non-cyclic boundaries
[707]179       IF ( .NOT. bc_lr_cyc )  THEN
[3182]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)
[1]182       ENDIF
[707]183       IF ( .NOT. bc_ns_cyc )  THEN
[3182]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,:)
[1]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
[4591]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)               )
[1]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
[4591]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)               )
[1]214             ENDDO
215          ENDDO
216       ENDDO
217
218!
219!--    Exchange of boundary values for p.
[667]220       CALL exchange_horiz( p, nbgp )
[1]221
222!
223!--    Boundary conditions top/bottom.
224!--    Bottom boundary
[4591]225       IF ( ibc_p_b == 1 )  THEN       ! Neumann
[1]226          p(nzb,:,:) = p(nzb+1,:,:)
[4591]227       ELSE                            ! Dirichlet
[1353]228          p(nzb,:,:) = 0.0_wp
[1]229       ENDIF
230
231!
232!--    Top boundary
[4591]233       IF ( ibc_p_t == 1 )  THEN       ! Neumann
[1]234          p(nzt+1,:,:) = p(nzt,:,:)
[4591]235       ELSE                            ! Dirichlet
[1353]236          p(nzt+1,:,:) = 0.0_wp
[1]237       ENDIF
238
239!
240!--    Horizontal (Neumann) boundary conditions in case of non-cyclic boundaries
[707]241       IF ( .NOT. bc_lr_cyc )  THEN
[3182]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)
[1]244       ENDIF
[707]245       IF ( .NOT. bc_ns_cyc )  THEN
[3182]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,:)
[1]248       ENDIF
249
[667]250
[1]251    ENDDO
252
253    DEALLOCATE( f1, f2, f3 )
254
[4457]255 END SUBROUTINE sor
Note: See TracBrowser for help on using the repository browser.