subroutine MATRIX_MATRIX &
   (ifile_in,ifile_out,name_in,name_out,NENS_START,NENS,MPIRANK,NPROC,jsta,N_LOC,X_ens)

!**********************************************************************
! *          
! * ABSTRACT: MATRIX * MATRIX product 
! *           (N_1 x NALL) * (NALL x NALL) = (N_1 x NALL)
! *               XIN      *  Matrix_NALL  =    XOUT
! *
! *  Use (N_1 x N_help) matrix = N_help columns of ouput matrix
! *
! *  INPUT: x_in        - columns of input matrix XIN=(N_1 x NALL)
! *                       [NALL columns of length (1:N_1)]
! *         X_ens       - Ensemble-space matrix (NALL x NALL)
! *
! * OUTPUT: x_out       - columns of output matrix XOUT=(N_1 x NALL)
! *                       [NALL columns of length (1:N_1)]
! *
! * 09/01/2003    ..... M. ZUPANSKI:
! * 10/17/2003    ..... M. ZUPANSKI: Generalized for arbitrary XIN, Matrix_eigen
! *
! **********************************************************************
!-----
     integer :: ifile_in,ifile_out
     integer :: NENS
     integer :: NENS_START
     integer :: NALL
     integer :: N_1,NN
     integer :: NPROC,MPIRANK
     integer :: N_LOC,jsta
     integer :: IO_max
     real,dimension(NENS_START:NENS,NENS_START:NENS) :: X_ens
     character*20 :: name_in,name_out
     character*30 :: filein,fileout

     real,allocatable::sum_Pa(:,:)
     real,allocatable::x_in(:)
     real,allocatable::x_out(:)
     integer,allocatable :: N_help(:)
!-----
! NAMELISTS
!-----
     NAMELIST /MEMORY/ MEM_PARAM
!=========================================================
        if(MPIRANK.eq.0) then
        write(*,*) "start MATRIX_MATRIX"
        end if

        NALL=NENS-NENS_START+1

!-- Parameter fitting the CPU size
!   (N_help*N_1) is largest memory space desired      
!   N_help_max = 1     => (N_LOC * NALL)=Num of I/Os, min memory (N_1)
!   N_help_max = N_LOC => (1 * NALL)=Num of I/Os, max memory (N_1 * N_LOC)
!   N_help_max         => N_LOC/N_help_max=Num of I/Os, memory (N_1 * N_help_max)

!- read memory size parameter
      rewind 14
      read(14,MEMORY) 

      N_help_max=min(MEM_PARAM,N_LOC)

      if(MOD(N_LOC,N_help_max).eq.0) then
       IO_max=N_LOC/N_help_max
       allocate(N_help(0:IO_max))
       N_help(0)=0
       do n=1,IO_max
       N_help(n)=N_help_max
       end do
      else
       IO_max=1+N_LOC/N_help_max
       allocate(N_help(0:IO_max))
       N_help(0)=0
       Nsum=0
       do n=1,IO_max-1
       N_help(n)=N_help_max
       Nsum=Nsum+N_help_max
       end do
       N_help(IO_max)=N_LOC-Nsum
      end if

      if(MPIRANK.eq.0) then
      write(*,*) MPIRANK," Number of I/O=",IO_max
      end if

  100    format(A,i4.4)

         write(filein,100) trim(name_in),jsta
      CLOSE(ifile_in)
      OPEN(UNIT=ifile_in,FILE=filein,FORM='UNFORMATTED',IOSTAT=IER)
      IF(IER.NE.0) WRITE(*,*) ifile_in,' OPEN UNIT ERROR IER=',IER
      READ(ifile_in) N_1
      CLOSE(ifile_in,STATUS='KEEP')

      allocate(x_in(1:N_1))
      allocate(x_out(1:N_1))

      do i=1,IO_max
      write(*,*) i," N_help=",N_help(i)
      end do
!==============start calculation===================
! jj=actual member
! ii=same member, but relative to the 1:N_help(n) loop
      n_read=0
      n_write=0

      j1=jsta
      do n=1,IO_max

      allocate(sum_Pa(1:N_1,1:N_help(n)))
      sum_Pa=0.0

      j1=j1+N_help(n-1)
      j2=j1+N_help(n)-1

         do i=NENS_START,NENS

         write(filein,100) trim(name_in),i
      CLOSE(ifile_in)
      OPEN(UNIT=ifile_in,FILE=filein,FORM='UNFORMATTED',IOSTAT=IER)
      IF(IER.NE.0) WRITE(*,*) ifile_in,' OPEN UNIT ERROR IER=',IER
      READ(ifile_in) NN
      if(NN.ne.N_1) then
      write(*,*) "MATRIX_MATRIX: PROBLEM WITH DIMENSIONS NN=",NN," N_1=",N_1
      end if
      READ(ifile_in) x_in
      n_read=n_read+1
      CLOSE(ifile_in,STATUS='KEEP')

      ii=0
       do jj=j1,j2
      ii=ii+1
       do m=1,N_1
       sum_Pa(m,ii)=sum_Pa(m,ii)+x_in(m)*X_ens(i,jj)
       end do
       end do

         end do

!---------------
        sumx11=0.0
        sumx22=0.0
        sumx12=0.0
        do is=1,N_1
        sumx11=sumx11+sum_Pa(is,1)*sum_Pa(is,1)
        sumx22=sumx22+sum_Pa(is,j2-j1+1)*sum_Pa(is,j2-j1+1)
        sumx12=sumx12+sum_Pa(is,1)*sum_Pa(is,j2-j1+1)
        end do
        write(*,*) j1,j2," MATRIX_MATRIX: sumx11,22,12=",sumx11,sumx22,sumx12
!---------------

      ii=0
       do jj=j1,j2
      ii=ii+1
      x_out(:)=sum_Pa(:,ii)
         write(fileout,100) trim(name_out),jj
      CLOSE(ifile_out)
      OPEN(UNIT=ifile_out,FILE=fileout,FORM='UNFORMATTED',IOSTAT=IER)
      IF(IER.NE.0) WRITE(*,*) ifile_out,' OPEN UNIT ERROR IER=',IER
      WRITE(ifile_out) N_1
      WRITE(ifile_out) x_out
      n_write=n_write+1
      CLOSE(ifile_out,STATUS='KEEP')

!--- statistics --------------------
                                                                                                  
        sumP=0.0
        do i=1,N_1
        sumP=sumP+x_out(i)**2
        end do
        sumP=sqrt(sumP/float(N_1))
      write(*,200) fileout,minval(x_out),maxval(x_out),sumP
!-----------------------------------

       end do       !jj

 200  format(a20,' min,max,avg=',3E15.5)

      deallocate(sum_Pa)

      end do        !n

      deallocate(x_in)
      deallocate(x_out)
      deallocate(N_help)

      write(*,*) "n_read=",n_read," n_write=",n_write

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

            return
            end subroutine MATRIX_MATRIX