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

 SUBROUTINE sor( d, ddzu, ddzw, p )

!------------------------------------------------------------------------------!
! Actual revisions:
! -----------------
! 
!
! Former revisions:
! -----------------
! $Id: sor.f90 4 2007-02-13 11:33:16Z raasch $
! 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: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) )

!
!-- 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, 0, 0 )

!
!--    Horizontal (Neumann) boundary conditions in case of non-cyclic boundaries
       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

!
!--    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, 0, 0 )

!
!--    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 ( 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

    DEALLOCATE( f1, f2, f3 )

 END SUBROUTINE sor