C The following program is used to test  if MPICH2's  predefined reduce operations support quadruple precision
      PROGRAM TEST

	IMPLICIT NONE

	INCLUDE 'MPIF.H'

C MPI function status return 
	INTEGER::ISTATUS(MPI_STATUS_SIZE),IERR
C My rank 
	INTEGER MYID

C Define the real precision
	INTEGER,PARAMETER::RP=SELECTED_REAL_KIND(33)  !one of 6,15,33
C MPI REAL DATA TYPE
      INTEGER MPI_REAL_RP

	REAL(KIND=RP) A,S

C Initilize the MPI enviroment
	CALL MPI_INIT(IERR)

      CALL MPI_TYPE_CREATE_F90_REAL(PRECISION(1.0_RP),MPI_UNDEFINED,
     &                              MPI_REAL_RP,IERR)

C Due the bug in mpich2, https://trac.mcs.anl.gov/projects/mpich2/ticket/1028
C use the following method to choose MPI REAL type
      IF(PRECISION(1.0_RP).EQ.6) THEN
	  MPI_REAL_RP=MPI_REAL
	ELSE IF(PRECISION(1.0_RP).EQ.15) THEN
	  MPI_REAL_RP=MPI_DOUBLE_PRECISION 
	ELSE IF(PRECISION(1.0_RP).EQ.33) THEN
	  MPI_REAL_RP=MPI_REAL16
	ELSE
	  WRITE(6,*)"MPI can not support specified single precision."
	  STOP
	END IF


C Get the rank ID
      CALL MPI_COMM_RANK(MPI_COMM_WORLD,MYID,IERR)

	A=MYID*1.0_RP

C Perform reduce operation
      CALL MPI_REDUCE(A,S,1,MPI_REAL_RP,MPI_SUM,0,MPI_COMM_WORLD,IERR)

      IF(MYID.EQ.0) THEN
	  PRINT*,'SUM=',S
	END IF

      
	CALL MPI_FINALIZE(IERR)


	END