program get_sqrtPx
!
!---   create ensemble perturbation (sqrtPa or sqrtPf) ----------
!
!   INPUT: x_cntrl = x_cntrl         
!          x       = ensemble forecast
!
!  OUTPUT: sqrtPx  = x - x_cntrl 
!
!--------------------------------------------------------------------

#ifdef MPI_USE
  include "mpif.h"
#endif

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

         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

     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 get_sqrtPx"
       end if

!-- read ensemble size

       REWIND 15
       READ(15,ENSEMBLE_SIZE)

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

!====================================================================
!-- 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')

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

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

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

!-- get the output string 

    call getenv("sqrtPx",Px)
          if(MPIRANK.eq.0) then
          write(*,*) "sqrtPx=",Px
          end if

     write(*,*) "         min,max x_cntrl=",minval(x_cntrl),maxval(x_cntrl)

      DO N=jsta,jend

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

     sqrtPx(:) = x(:) - x_cntrl(:)

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

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

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

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

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

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

#ifdef MPI_USE
      CALL MPI_FINALIZE(IERR)
#endif


stop
end program get_sqrtPx