SUBROUTINE sor( d, ddzu, ddzw, p )

!------------------------------------------------------------------------------!
! Current revisions:
! -----------------
! 
!
! Former revisions:
! -----------------
! $Id: sor.f90 708 2011-03-29 12:34:54Z gryschka $
!
! 707 2011-03-29 11:39:40Z raasch
! bc_lr/ns replaced by bc_lr/ns_cyc
!
! 667 2010-12-23 12:06:00Z suehring/gryschka
! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.
! Call of exchange_horiz are modified.
! bug removed in declaration of ddzw(), nz replaced by nzt+1
!
! 75 2007-03-22 09:54:05Z raasch
! 2nd+3rd argument removed from exchange horiz
!
! RCS Log replace by Id keyword, revision history cleaned up
!
! 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.1  1997/08/11 06:25:56  raasch
! Initial revision
!
!
! Description:
! ------------
! Solve the Poisson-equation with the SOR-Red/Black-scheme.
!------------------------------------------------------------------------------!

    USE grid_variables
    USE indices
    USE pegrid
    USE control_parameters

    IMPLICIT NONE

    INTEGER ::  i, j, k, n, nxl1, nxl2, nys1, nys2
    REAL    ::  ddzu(1:nz+1), ddzw(1:nzt+1)
    REAL    ::  d(nzb+1:nzt,nys:nyn,nxl:nxr),         &
                p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)
    REAL, DIMENSION(:), ALLOCATABLE ::  f1, f2, f3

    ALLOCATE( f1(1:nz), f2(1:nz), f3(1:nz) )

!
!-- Compute pre-factors.
    DO  k = 1, nz
         f2(k) = ddzu(k+1) * ddzw(k)
         f3(k) = ddzu(k)   * ddzw(k)
         f1(k) = 2.0 * ( ddx2 + ddy2 ) + f2(k) + f3(k)
    ENDDO

!
!-- Limits for RED- and BLACK-part.
    IF ( MOD( nxl , 2 ) == 0 )  THEN
       nxl1 = nxl
       nxl2 = nxl + 1
    ELSE
       nxl1 = nxl + 1
       nxl2 = nxl
    ENDIF
    IF ( MOD( nys , 2 ) == 0 )  THEN
       nys1 = nys
       nys2 = nys + 1
    ELSE
       nys1 = nys + 1
       nys2 = nys
    ENDIF

    DO  n = 1, n_sor

!
!--    RED-part
       DO  i = nxl1, nxr, 2
          DO  j = nys2, nyn, 2
             DO  k = nzb+1, nzt
                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
       ENDDO

       DO  i = nxl2, nxr, 2
          DO  j = nys1, nyn, 2
             DO  k = nzb+1, nzt
                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
       ENDDO

!
!--    Exchange of boundary values for p.
       CALL exchange_horiz( p, nbgp )

!
!--    Horizontal (Neumann) boundary conditions in case of non-cyclic boundaries
       IF ( .NOT. bc_lr_cyc )  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 ( .NOT. bc_ns_cyc )  THEN
          IF ( inflow_n .OR. outflow_n )  p(:,nyn+1,:) = p(:,nyn,:)
          IF ( inflow_s .OR. outflow_s )  p(:,nys-1,:) = p(:,nys,:)
       ENDIF

!
!--    BLACK-part
       DO  i = nxl1, nxr, 2
          DO  j = nys1, nyn, 2
             DO  k = nzb+1, nzt
                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
       ENDDO

       DO  i = nxl2, nxr, 2
          DO  j = nys2, nyn, 2
             DO  k = nzb+1, nzt
                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
       ENDDO

!
!--    Exchange of boundary values for p.
       CALL exchange_horiz( p, nbgp )

!
!--    Boundary conditions top/bottom.
!--    Bottom boundary
       IF ( ibc_p_b == 1 )  THEN       !       Neumann
          p(nzb,:,:) = p(nzb+1,:,:)
       ELSE                            !       Dirichlet
          p(nzb,:,:) = 0.0
       ENDIF

!
!--    Top boundary
       IF ( ibc_p_t == 1 )  THEN                 !  Neumann
          p(nzt+1,:,:) = p(nzt,:,:)
       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 ( inflow_l .OR. outflow_l )  p(:,:,nxl-1) = p(:,:,nxl)
          IF ( inflow_r .OR. outflow_r )  p(:,:,nxr+1) = p(:,:,nxr)
       ENDIF
       IF ( .NOT. bc_ns_cyc )  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

    DEALLOCATE( f1, f2, f3 )

 END SUBROUTINE sor