C======================================================================
C NPAC $Id: find_rhs_psi.f90,v 1.4 1995/11/21 20:42:23 hwyau Exp $
C======================================================================
C
      SUBROUTINE find_rhs_psi(alpha,h2,psi,zeta,rhs,N)
C _____________________________________________________________________
C|                                                                     |
C| Subroutine, to find the RHS for ADI to solve the Poisson equation   |
C| with time-dependent component.                                      |
C|                                                                     |
C|    RHS = A*psi - h2*zeta                                            |
C| Where the matrix A is tridiagonal, with unit off diagonal elements. |
C|                                                                     |
C| H W Yau.  17th of April, 1995.                                      |
C| Northeast Parallel Architectures Center.                            |
C| Syracuse University.                                                |
C|_____________________________________________________________________|
C|                                                                     |
C| Edit record:                                                        |
C| 17/Apr/1995: First edit.                                            |
C| 19/Apr/1995: Stuck in code.                                         |
C| 15/Nov/1995: Rewrite, to avoid INDEPENDENT.                         |
C|              Added HPF directives.                                  |
C| 17/Nov/1995: Option to not update boundary elements.                |
C| 21/Nov/1995: Correction to handling of boundary elements.           |
C|_____________________________________________________________________|
C
      IMPLICIT NONE
C
C .. Arguments ..
      INTEGER,INTENT(IN)              :: N
      REAL,INTENT(IN)                 :: alpha,h2
      REAL,DIMENSION(:,:),INTENT(IN)  :: psi,zeta
      REAL,DIMENSION(:,:),INTENT(OUT) :: rhs
C
C .. HPF directives ..
!HPF$ DISTRIBUTE *(*,BLOCK) :: psi,zeta,rhs
C
C .. Local variables ..
      INTEGER i,k,j1,jN
      REAL b
C______________________________________________________________________
C
C Begin Executable
C______________________________________________________________________
C
      j1 = 2                        ! Start element to update
      jN = N-1                      ! End index to update
      b = alpha - 2.0               ! Diagonal elements
      DO k=j1,jN
        rhs(k,j1:jN) = psi(k-1,j1:jN) + b*psi(k,j1:jN) + psi(k+1,j1:jN)
     $               - h2*zeta(k,j1:jN)
      END DO
C
C Set boundary conditions
      rhs(1,:) = psi(1,:)
      rhs(N,:) = psi(N,:)
      DO k=j1,jN
        rhs(k,1) = psi(k,1)
        rhs(k,N) = psi(k,N)
      END DO
C
      RETURN
C
      END SUBROUTINE find_rhs_psi
C
C======================================================================

