program zcntrl_mean_categ
!
!---   calculate the ensemble mean ---------------------
!
!   INPUT: x       = x                           !!! ensemble members
!  OUTPUT: x_mean  = x_mean                      !!! ensemble mean
!
!--------------------------------------------------------------------

#ifdef MPI_USE
 include "mpif.h"
#endif

   integer,parameter::iptrb=20            ! input ensemble member 
   integer,parameter::imean=51            ! output ensemble mean (x1d_cntrl)
!-----
   character*12 filename
         real,dimension(:),allocatable :: x_mean,x
         real,dimension(:),allocatable :: sum_loc,sum_gbl
         integer,allocatable :: jlen(:),jdisp(:)

     integer :: N_LOC,jsta,jend
     integer :: NPROC,MPIRANK
     integer :: NENS
     integer :: NENS_START
     integer :: NALL
     integer :: ncateg,N_start,N_end

     NAMELIST /ENSEMBLE_SIZE/ NENS_START,NENS
!====================================================================

#ifdef MPI_USE
!
!     START MPI
!
      CALL MPI_INIT(IERR)
      IF (IERR .NE. MPI_SUCCESS) STOP 'FAILED TO INIT MPI'
      CALL MPI_COMM_SIZE(MPI_COMM_WORLD, NPROC, IERR)
      CALL MPI_COMM_RANK(MPI_COMM_WORLD, MPIRANK, IERR)
      CALL MPI_BARRIER(MPI_COMM_WORLD,IERR)
#else
      write(*,*) "This is a NO MPI run"
      MPIRANK=0
      NPROC=1
#endif

       if(MPIRANK.eq.0) then
       write(*,*) "start mean"
       end if

!-- read ensemble size

       REWIND 15
       READ(15,ENSEMBLE_SIZE)


!-- this is the mean
        NALL=NENS-NENS_START+1

          if(MPIRANK.eq.0) then
          write(*,*) "mean: number of ensemble members=",NALL
          end if

!--Read number of categories of matrix A-------
                                                                                
      read(17,*) ncateg
                                                                                
      N_start=NENS_START
      N_end=NALL*ncateg
!====================================================================
!==
!==  define local dimensions and indexes
!==
       allocate(jdisp(0:NPROC-1)) ; allocate(jlen(0:NPROC-1))

            do irank=0,NPROC-1
      CALL para_range(NENS_START,NENS,NPROC,irank,jsta,jend)
            jlen(irank)=jend-jsta+1
            jdisp(irank)=jsta-NENS_START
            end do
#ifdef MPI_USE
      CALL MPI_BARRIER(MPI_COMM_WORLD,IERR)
#endif
      CALL para_range(NENS_START,NENS,NPROC,MPIRANK,jsta,jend)

!==
         N_LOC=jlen(MPIRANK)
!==
!
         if(MPIRANK.eq.0) then
           do nn=0,NPROC-1
           write(*,*) "nn=",nn," jdisp=",jdisp(nn)," jlen=",jlen(nn)
           end do
          end if
!=================================================================
!   idim=NALL
   idim=N_end-N_start+1
   allocate(x(1:idim))
   allocate(sum_loc(1:idim))
     sum_loc(:)=0.0

      DO N=jsta,jend

      WRITE(filename,1000)  N
      CLOSE(iptrb)
      OPEN(iptrb,file=filename,status='unknown',form='unformatted')
      READ(iptrb) x
      CLOSE(iptrb,status='KEEP')

      do iloc=1,idim
     sum_loc(iloc)=sum_loc(iloc) + x(iloc)
      end do

      END DO

 1000 format('z_cntrl.',i4.4)
      write(*,*) N," sum_loc=",minval(sum_loc),maxval(sum_loc)
!--------------------------------------------
!!! get allredauce sum: x_mean vs. sum_loc 

   allocate(sum_gbl(1:idim))
   sum_gbl=0.0

#ifdef MPI_USE
         CALL MPI_REDUCE &
 (sum_loc,sum_gbl,idim,MPI_REAL,MPI_SUM,0,MPI_COMM_WORLD,IERR)
#else
  sum_gbl(:)=sum_loc(:)
#endif

  deallocate(sum_loc)
  deallocate(x)

!--------------------------------------------

       if(MPIRANK.eq.0) then

   allocate(x_mean(1:idim))

      x_mean(:)=sum_gbl(:)/float(NALL)

      CLOSE(imean)
      OPEN(imean,file='z_cntrl',status='unknown',form='unformatted')
      WRITE(imean) x_mean
      CLOSE(imean,status='KEEP')

      deallocate(x_mean)
       write(*,*) "end mean"
       end if

   deallocate(sum_gbl)

!====================================================================
#ifdef MPI_USE
      CALL MPI_FINALIZE(IERR)
#endif


stop
end program zcntrl_mean_categ