C======================================================================
C NPAC $Id: main.f90,v 1.2 1996/05/30 17:49:30 hwyau Exp $
C======================================================================
C
      PROGRAM adi
C _____________________________________________________________________
C|                                                                     |
C| Driver for 2 dimensional ADI program, in HPF.                       |
C| Thus far, solves a simple Poisson equation plus a time-dependent    |
C| diffusive term:                                                     |
C|                                                                     |
C| \frac{\partial\psi}{\partial t} =  \nabla^2\psi - \zeta             |
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: Made into a real program.                              |
C| 15/Nov/1995: Added statements to print the results.                 |
C|              Added TRANSPOSE() function.                            |
C|              Made use of parallelised subroutines, to replace use   |
C|              of INDEPENDENT.                                        |
C|              Added HPF directives.                                  |
C| 21/Nov/1995: Works, for Laplace equation.                           |
C| 29/May/1996: Increased problem size to obtain timings.              |
C| 30/May/1996: REU version.                                           |
C|_____________________________________________________________________|
C
      IMPLICIT NONE
C
C The size of the problem is NN*NN
      INTEGER, PARAMETER :: NN = 256
C
C Mesh and time step parameters.
      REAL, PARAMETER    :: DT = 0.1
      REAL, PARAMETER    :: H  = 0.2
C
C Standard out.
      INTEGER, PARAMETER :: STDOUT = 6
C
C .. Define interface for subroutines ..
      INTERFACE
        SUBROUTINE tridiag_psi(alpha,beta,rhs,u,temp,N)
          INTEGER,INTENT(IN)              :: N
          REAL,INTENT(IN)                 :: alpha
          REAL,DIMENSION(:,:),INTENT(IN)  :: rhs
          REAL,DIMENSION(:,:),INTENT(OUT) :: u,temp,beta
!HPF$     DISTRIBUTE *(*,BLOCK) :: beta,rhs,u,temp
        END SUBROUTINE tridiag_psi
C
        SUBROUTINE find_rhs_psi(alpha,h2,psi,zeta,rhs,N)
          INTEGER,INTENT(IN)              :: N
          REAL,INTENT(IN)                 :: alpha,h2
          REAL,DIMENSION(:,:),INTENT(IN)  :: psi,zeta
          REAL,DIMENSION(:,:),INTENT(OUT) :: rhs
!HPF$     DISTRIBUTE *(*,BLOCK) :: psi,zeta,rhs
        END SUBROUTINE find_rhs_psi
      END INTERFACE
C
C .. Arrays and scalars ..
      INTEGER               :: i,j
      INTEGER               :: iTime
C
C The number of time steps to iterate for.
      INTEGER               :: NTime
      REAL                  :: alpha,h2
      REAL,DIMENSION(NN,NN) :: beta,psi,zeta,rhs,wspace
C
C .. HPF dirctives to distribute arrays ..
!HPF$ PROCESSORS             :: procs(NUMBER_OF_PROCESSORS())
!HPF$ TEMPLATE,DIMENSION(NN) :: tplate
!HPF$ DISTRIBUTE (BLOCK) ONTO procs :: tplate
!HPF$ ALIGN beta(*,:)   WITH tplate(:)
!HPF$ ALIGN psi(*,:)    WITH tplate(:)
!HPF$ ALIGN zeta(*,:)   WITH tplate(:)
!HPF$ ALIGN rhs(*,:)    WITH tplate(:)
!HPF$ ALIGN wspace(*,:) WITH tplate(:)
C
C______________________________________________________________________
C
C Begin Executable
C______________________________________________________________________
C
C .. Initialise scalars ..
      h2    = H*H
      alpha = 2.0*h2/DT      
      NTime = 10
C
C .. Initialise arrays ..
      rhs    = 1.0
      wspace = 0.0
      psi    = 5.0
      psi(2:NN-1,2:NN-1) = 2.0
C
C .. Reduce to Laplace equation.
      zeta = 0.0
C
C .. Solve, over a number of time steps ..
      DO iTime = 1,NTime
C
C .. Solve in x-direction.
        CALL tridiag_psi(alpha,beta,rhs,psi,wspace,NN)
        CALL find_rhs_psi(alpha,h2,psi,zeta,rhs,NN)
C
        psi    = TRANSPOSE(psi)
        zeta   = TRANSPOSE(zeta)
        rhs    = TRANSPOSE(rhs)
C
C .. Solve in y-direction.
        CALL tridiag_psi(alpha,beta,rhs,psi,wspace,NN)
        CALL find_rhs_psi(alpha,h2,psi,zeta,rhs,NN)
C
        psi    = TRANSPOSE(psi)
        zeta   = TRANSPOSE(zeta)
        rhs    = TRANSPOSE(rhs)
C
        WRITE(UNIT=STDOUT,FMT=1001) iTime
CC        DO i=1,NN
CC          WRITE(UNIT=STDOUT,FMT=1002) (psi(i,j),j=1,NN)
CC        END DO
      END DO ! iTime
C
CC      WRITE(UNIT=STDOUT,FMT=*) 'Final Result: for psi()'
CC      DO i=1,NN
CC        WRITE(UNIT=STDOUT,FMT=1002) (psi(i,j),j=1,NN)
CC      END DO
C
      STOP
C
 1001 FORMAT('Status: iTime=',I6)
 1002 FORMAT(12F9.3)
C
      END
C
C======================================================================
