program get_x1d
!
!---   create ensemble state (x1d_...) ----------
!
!   INPUT: x_cntrl = x_cntrl         
!          sqrtPx  = sqrtPa or sqrtPf
!
!  OUTPUT: x       = x_cntrl + sqrtPx - ensemble state
!
!--------------------------------------------------------------------

#ifdef XYZ
  include "mpif.h"
#endif

   integer,parameter::icntrl=20            ! cntrol 
   integer,parameter::isqrtPx=21           ! perturbation
   integer,parameter::iptrb=50             ! ensemble state

         integer idim
!-----
         real,dimension(:),allocatable :: x_cntrl
         real,dimension(:),allocatable :: sqrtPx,x
         integer,allocatable :: jlen(:),jdisp(:)

     integer :: N_LOC,jsta,jend
     integer :: NPROC,MPIRANK
     integer :: NENS
     integer :: NENS_START
     integer :: NALL
     character*8 filename
     character*11 Pxname
     character*7 Px
     character*4 stattype

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


#ifdef XYZ
!
!     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 get_x1d"
       end if

!-- read ensemble size

       REWIND 15
       READ(15,ENSEMBLE_SIZE)

          if(MPIRANK.eq.0) then
          write(*,*) "get_x1d: NENS_START, NENS=",NENS_START,NENS
          end if

!- Total number of degrees of freedom
        NALL=NENS-NENS_START+1
!====================================================================
!-- read control variable

      CLOSE(icntrl)
      OPEN(icntrl,file='x1d_cntrl',status='unknown',form='unformatted')
      READ(icntrl) idim
      allocate(x_cntrl(1:idim))
      READ(icntrl) x_cntrl
      CLOSE(icntrl,status='KEEP')

!====================================================================
!==
!==  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 XYZ
      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
!====================================================================
   allocate(sqrtPx(1:idim))
   allocate(x(1:idim))

 1000 format('x1d_',i4.4)
 2000 format(a7,i4.4)

!-- get the output string 

    call getenv("sqrtP",Px)

      DO N=jsta,jend

!-- read perturbation (sqrtPx)
      WRITE(Pxname,2000)  Px,N
      CLOSE(isqrtPx)
      OPEN(isqrtPx,file=Pxname,status='unknown',form='unformatted')
      READ(isqrtPx) idim
      READ(isqrtPx) sqrtPx
      CLOSE(isqrtPx,status='KEEP')

      write(*,*) N," sqrtPx min,max=",minval(sqrtPx),maxval(sqrtPx)
      write(*,*) N," x_cntrl min,max=",minval(x_cntrl),maxval(x_cntrl)
      x(:) = x_cntrl(:) + sqrtPx(:)
      write(*,*) N," x min,max=",minval(x),maxval(x)

!-- write ensemble state 
      WRITE(filename,1000)  N
      CLOSE(iptrb)
      OPEN(iptrb,file=filename,status='unknown',form='unformatted')
      WRITE(iptrb) idim
      WRITE(iptrb) x
      CLOSE(iptrb,status='KEEP')

      END DO
!--------------------------------------------

      deallocate(x)
      deallocate(x_cntrl)
      deallocate(sqrtPx)

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

!====================================================================

#ifdef XYZ
      CALL MPI_FINALIZE(IERR)
#endif


stop
end program get_x1d