Index: palm/trunk/SOURCE/sor.f90
===================================================================
--- palm/trunk/SOURCE/sor.f90	(revision 4901)
+++ palm/trunk/SOURCE/sor.f90	(revision 1)
@@ -1,113 +1,59 @@
-!> @file sor.f90
-!--------------------------------------------------------------------------------------------------!
-! This file is part of the PALM model system.
-!
-! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
-! Public License as published by the Free Software Foundation, either version 3 of the License, or
-! (at your option) any later version.
-!
-! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
-! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
-! Public License for more details.
-!
-! You should have received a copy of the GNU General Public License along with PALM. If not, see
-! <http://www.gnu.org/licenses/>.
-!
-! Copyright 1997-2021 Leibniz Universitaet Hannover
-!--------------------------------------------------------------------------------------------------!
-!
-!
-! Current revisions:
+ SUBROUTINE sor( d, ddzu, ddzw, p )
+
+!------------------------------------------------------------------------------!
+! Actual revisions:
 ! -----------------
-!
+! 
 !
 ! Former revisions:
 ! -----------------
-! $Id$
-! File re-formatted to follow the PALM coding standard
-!
-!
-! 4457 2020-03-11 14:20:43Z raasch
-! Use statement for exchange horiz added
-!
-! 4360 2020-01-07 11:25:50Z suehring
-! Corrected "Former revisions" section
-!
-! 3655 2019-01-07 16:51:22Z knoop
-! Rename variables in mesoscale-offline nesting mode
+! $Log: sor.f90,v $
+! Revision 1.9  2005/03/26 21:02:23  raasch
+! Implementation of non-cyclic (Neumann) horizontal boundary conditions,
+! dx2,dy2 replaced by ddx2,ddy2
+!
+! Revision 1.8  2001/03/30 07:52:44  raasch
+! Arrays f1-f3 changed from dummy arguments to allocatable arrays,
+! Translation of remaining German identifiers (variables, subroutines, etc.)
+!
+! Revision 1.7  2001/01/30 21:41:25  letzel
+! All comments translated into English.
+!
+! Revision 1.6  2001/01/22 08:06:14  raasch
+! Module test_variables removed
+!
+! Revision 1.5  1998/07/06 12:33:13  raasch
+! + USE test_variables
+!
+! Revision 1.4  1997/09/12 06:30:12  raasch
+! Randbedingungen umbenannt
+!
+! Revision 1.3  1997/09/09 08:30:44  raasch
+! Kehrwerte der Gitterweiten implementiert
+!
+! Revision 1.2  1997/08/29 09:00:58  raasch
+! omega --> omega_sor
 !
 ! Revision 1.1  1997/08/11 06:25:56  raasch
 ! Initial revision
 !
-!--------------------------------------------------------------------------------------------------!
+!
 ! Description:
 ! ------------
-!> Solve the Poisson-equation with the SOR-Red/Black-scheme.
-!--------------------------------------------------------------------------------------------------!
- SUBROUTINE sor( d, ddzu, ddzw, p )
-
-    USE arrays_3d,                                                                                 &
-        ONLY:  rho_air,                                                                            &
-               rho_air_zw
-
-    USE control_parameters,                                                                        &
-        ONLY:  bc_dirichlet_l,                                                                     &
-               bc_dirichlet_n,                                                                     &
-               bc_dirichlet_r,                                                                     &
-               bc_dirichlet_s,                                                                     &
-               bc_lr_cyc,                                                                          &
-               bc_ns_cyc,                                                                          &
-               bc_radiation_l,                                                                     &
-               bc_radiation_n,                                                                     &
-               bc_radiation_r,                                                                     &
-               bc_radiation_s,                                                                     &
-               ibc_p_b,                                                                            &
-               ibc_p_t,                                                                            &
-               n_sor,                                                                              &
-               omega_sor
-
-    USE exchange_horiz_mod,                                                                        &
-        ONLY:  exchange_horiz
-
-    USE grid_variables,                                                                            &
-        ONLY:  ddx2,                                                                               &
-               ddy2
-
-    USE indices,                                                                                   &
-        ONLY:  nbgp,                                                                               &
-               nxl,                                                                                &
-               nxlg,                                                                               &
-               nxr,                                                                                &
-               nxrg,                                                                               &
-               nyn,                                                                                &
-               nyng,                                                                               &
-               nys,                                                                                &
-               nysg,                                                                               &
-               nz,                                                                                 &
-               nzb,                                                                                &
-               nzt
-
-    USE kinds
+! Solve the Poisson-equation with the SOR-Red/Black-scheme.
+!-------------------------------------------------------------------------------!
+
+    USE grid_variables
+    USE indices
+    USE pegrid
+    USE control_parameters
 
     IMPLICIT NONE
 
-    INTEGER(iwp) ::  i     !<
-    INTEGER(iwp) ::  j     !<
-    INTEGER(iwp) ::  k     !<
-    INTEGER(iwp) ::  n     !<
-    INTEGER(iwp) ::  nxl1  !<
-    INTEGER(iwp) ::  nxl2  !<
-    INTEGER(iwp) ::  nys1  !<
-    INTEGER(iwp) ::  nys2  !<
-
-    REAL(wp) ::  ddzu(1:nz+1)   !<
-    REAL(wp) ::  ddzw(1:nzt+1)  !<
-
-    REAL(wp) ::  d(nzb+1:nzt,nys:nyn,nxl:nxr)      !<
-    REAL(wp) ::  p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  !<
-
-    REAL(wp), DIMENSION(:), ALLOCATABLE ::  f1  !<
-    REAL(wp), DIMENSION(:), ALLOCATABLE ::  f2  !<
-    REAL(wp), DIMENSION(:), ALLOCATABLE ::  f3  !<
+    INTEGER ::  i, j, k, n, nxl1, nxl2, nys1, nys2
+    REAL    ::  ddzu(1:nz+1), ddzw(1:nz)
+    REAL    ::  d(nzb+1:nzt,nys:nyn,nxl:nxr),         &
+                p(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
+    REAL, DIMENSION(:), ALLOCATABLE ::  f1, f2, f3
 
     ALLOCATE( f1(1:nz), f2(1:nz), f3(1:nz) )
@@ -116,7 +62,7 @@
 !-- Compute pre-factors.
     DO  k = 1, nz
-         f2(k) = ddzu(k+1) * ddzw(k) * rho_air_zw(k)
-         f3(k) = ddzu(k)   * ddzw(k) * rho_air_zw(k-1)
-         f1(k) = 2.0_wp * ( ddx2 + ddy2 ) * rho_air(k) + f2(k) + f3(k)
+         f2(k) = ddzu(k+1) * ddzw(k)
+         f3(k) = ddzu(k)   * ddzw(k)
+         f1(k) = 2.0 * ( ddx2 + ddy2 ) + f2(k) + f3(k)
     ENDDO
 
@@ -145,11 +91,11 @@
           DO  j = nys2, nyn, 2
              DO  k = nzb+1, nzt
-                p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * (                                        &
-                           rho_air(k) * ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +                       &
-                           rho_air(k) * ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +                       &
-                           f2(k) * p(k+1,j,i)                              +                       &
-                           f3(k) * p(k-1,j,i)                              -                       &
-                           d(k,j,i)                                        -                       &
-                           f1(k) * p(k,j,i)               )
+                p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * (            &
+                               ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +    &
+                               ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +    &
+                               f2(k) * p(k+1,j,i)                 +    &
+                               f3(k) * p(k-1,j,i)                 -    &
+                               d(k,j,i)                           -    &
+                               f1(k) * p(k,j,i)           )
              ENDDO
           ENDDO
@@ -159,11 +105,11 @@
           DO  j = nys1, nyn, 2
              DO  k = nzb+1, nzt
-                p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * (                                        &
-                           rho_air(k) * ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +                       &
-                           rho_air(k) * ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +                       &
-                           f2(k) * p(k+1,j,i)                              +                       &
-                           f3(k) * p(k-1,j,i)                              -                       &
-                           d(k,j,i)                                        -                       &
-                           f1(k) * p(k,j,i)               )
+                p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * (            &
+                               ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +    &
+                               ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +    &
+                               f2(k) * p(k+1,j,i)                 +    &
+                               f3(k) * p(k-1,j,i)                 -    &
+                               d(k,j,i)                           -    &
+                               f1(k) * p(k,j,i)           )
              ENDDO
           ENDDO
@@ -172,16 +118,15 @@
 !
 !--    Exchange of boundary values for p.
-       CALL exchange_horiz( p, nbgp )
-
+       CALL exchange_horiz( p, 0, 0 )
 
 !
 !--    Horizontal (Neumann) boundary conditions in case of non-cyclic boundaries
-       IF ( .NOT. bc_lr_cyc )  THEN
-          IF ( bc_dirichlet_l  .OR.  bc_radiation_l )  p(:,:,nxl-1) = p(:,:,nxl)
-          IF ( bc_dirichlet_r  .OR.  bc_radiation_r )  p(:,:,nxr+1) = p(:,:,nxr)
-       ENDIF
-       IF ( .NOT. bc_ns_cyc )  THEN
-          IF ( bc_dirichlet_n  .OR.  bc_radiation_n )  p(:,nyn+1,:) = p(:,nyn,:)
-          IF ( bc_dirichlet_s  .OR.  bc_radiation_s )  p(:,nys-1,:) = p(:,nys,:)
+       IF ( bc_lr /= 'cyclic' )  THEN
+          IF ( inflow_l .OR. outflow_l )  p(:,:,nxl-1) = p(:,:,nxl)
+          IF ( inflow_r .OR. outflow_r )  p(:,:,nxr+1) = p(:,:,nxr)
+       ENDIF
+       IF ( bc_ns /= 'cyclic' )  THEN
+          IF ( inflow_n .OR. outflow_n )  p(:,nyn+1,:) = p(:,nyn,:)
+          IF ( inflow_s .OR. outflow_s )  p(:,nys-1,:) = p(:,nys,:)
        ENDIF
 
@@ -191,11 +136,11 @@
           DO  j = nys1, nyn, 2
              DO  k = nzb+1, nzt
-                p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * (                                        &
-                           rho_air(k) * ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +                       &
-                           rho_air(k) * ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +                       &
-                           f2(k) * p(k+1,j,i)                              +                       &
-                           f3(k) * p(k-1,j,i)                              -                       &
-                           d(k,j,i)                                        -                       &
-                           f1(k) * p(k,j,i)               )
+                p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * (            &
+                               ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +    &
+                               ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +    &
+                               f2(k) * p(k+1,j,i)                 +    &
+                               f3(k) * p(k-1,j,i)                 -    &
+                               d(k,j,i)                           -    &
+                               f1(k) * p(k,j,i)           )
              ENDDO
           ENDDO
@@ -205,11 +150,11 @@
           DO  j = nys2, nyn, 2
              DO  k = nzb+1, nzt
-                p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * (                                        &
-                           rho_air(k) * ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +                       &
-                           rho_air(k) * ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +                       &
-                           f2(k) * p(k+1,j,i)                              +                       &
-                           f3(k) * p(k-1,j,i)                              -                       &
-                           d(k,j,i)                                        -                       &
-                           f1(k) * p(k,j,i)               )
+                p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * (            &
+                               ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +    &
+                               ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +    &
+                               f2(k) * p(k+1,j,i)                 +    &
+                               f3(k) * p(k-1,j,i)                 -    &
+                               d(k,j,i)                           -    &
+                               f1(k) * p(k,j,i)           )
              ENDDO
           ENDDO
@@ -218,34 +163,41 @@
 !
 !--    Exchange of boundary values for p.
-       CALL exchange_horiz( p, nbgp )
+       CALL exchange_horiz( p, 0, 0 )
 
 !
 !--    Boundary conditions top/bottom.
 !--    Bottom boundary
-       IF ( ibc_p_b == 1 )  THEN       ! Neumann
+       IF ( ibc_p_b == 1 )  THEN
+!
+!--       Neumann
           p(nzb,:,:) = p(nzb+1,:,:)
-       ELSE                            ! Dirichlet
-          p(nzb,:,:) = 0.0_wp
+       ELSE
+!
+!--       Dirichlet
+          p(nzb,:,:) = 0.0
        ENDIF
 
 !
 !--    Top boundary
-       IF ( ibc_p_t == 1 )  THEN       ! Neumann
+       IF ( ibc_p_t == 1 )  THEN
+!
+!--       Neumann
           p(nzt+1,:,:) = p(nzt,:,:)
-       ELSE                            ! Dirichlet
-          p(nzt+1,:,:) = 0.0_wp
+       ELSE
+!
+!--       Dirichlet
+          p(nzt+1,:,:) = 0.0
        ENDIF
 
 !
 !--    Horizontal (Neumann) boundary conditions in case of non-cyclic boundaries
-       IF ( .NOT. bc_lr_cyc )  THEN
-          IF ( bc_dirichlet_l  .OR.  bc_radiation_l )  p(:,:,nxl-1) = p(:,:,nxl)
-          IF ( bc_dirichlet_r  .OR.  bc_radiation_r )  p(:,:,nxr+1) = p(:,:,nxr)
-       ENDIF
-       IF ( .NOT. bc_ns_cyc )  THEN
-          IF ( bc_dirichlet_n  .OR.  bc_radiation_n )  p(:,nyn+1,:) = p(:,nyn,:)
-          IF ( bc_dirichlet_s  .OR.  bc_radiation_s )  p(:,nys-1,:) = p(:,nys,:)
-       ENDIF
-
+       IF ( bc_lr /= 'cyclic' )  THEN
+          IF ( inflow_l .OR. outflow_l )  p(:,:,nxl-1) = p(:,:,nxl)
+          IF ( inflow_r .OR. outflow_r )  p(:,:,nxr+1) = p(:,:,nxr)
+       ENDIF
+       IF ( bc_ns /= 'cyclic' )  THEN
+          IF ( inflow_n .OR. outflow_n )  p(:,nyn+1,:) = p(:,nyn,:)
+          IF ( inflow_s .OR. outflow_s )  p(:,nys-1,:) = p(:,nys,:)
+       ENDIF
 
     ENDDO
