program transform
!**********************************************************************
! *              .      .    .
! * PROGRAM:    TRANSFORM
! * PRGMMR:     M. ZUPANSKI     ORG: CIRA/CSU   DATE: 2003-09-09
! *
! * ABSTRACT:   TRANSFORMATION OF CONTROL VARIABLE FROM ENSEMBLE SPACE TO MODEL SPACE (1-d)
! *
! * PROGRAM LOG:
! *
! * 09/09/2003    ..... M. ZUPANSKI:
! *
! **********************************************************************
!----------------------------------
!  INPUT: 
!         sqrt(Pf) - fcst_xxxx
!         sqrt(I+A)**(-1/2) - sqrtIPA_inv
!         z - control variable in ensemble space
!
!  OUTPUT:
!         x - 1-d control variable in model space
!
!----------------------------------

#ifdef MPI_USE
 include "mpif.h"
#endif
!-----

     integer,parameter::xdim=20              ! dimensions (x1d_cntrl)
     integer,parameter::Pa_file=22        ! ensemble fcst ptrb file #
     integer,parameter::icntrl=23         ! control variable in ensemble space
     integer,parameter::imodelb=24        ! first guess control variable in model space
     integer,parameter::imodel=51         ! control variable in model space
!-----
     integer :: NENS
     integer :: NENS_START
     integer :: NALL
     integer :: N_model
     integer :: N_LOC,jsta,jend
     integer :: NPROC,MPIRANK
     character*20 :: namein
     character*30 :: filein

     integer,allocatable :: jlen(:),jdisp(:)

     real,dimension(:),allocatable::z
     real,dimension(:),allocatable::sum_loc,sumx
     real,dimension(:),allocatable::x,x_b
!-----
!
!      DECLARE NAMELIST
!
     NAMELIST /ENSEMBLE_SIZE/ NENS_START,NENS

!==============start calculation===================

#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
!================================================

       REWIND 15
       READ(15,ENSEMBLE_SIZE)

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

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

!-- input filename
      write(namein,101)
 101  format('sqrtPa_')

!-- model dimensions (N_model)
      CLOSE(xdim)
      OPEN(xdim,file='x1d_cntrl',status='unknown',form='unformatted')
      READ(xdim) N_model
      CLOSE(xdim,status='KEEP')
       if(MPIRANK.eq.0) then
       write(*,*) "transform: N_model=",N_model
       end if

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

!==
!==  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)
!==
!
           write(*,*) "NALL,NPROC,MPIRANK,jsta,jend=",NALL,NPROC,MPIRANK,jsta,jend
         if(MPIRANK.eq.0) then
           do nn=0,NPROC-1
           write(*,*) "nn=",nn," jdisp=",jdisp(nn)," jlen=",jlen(nn)
           end do
          end if
!====================================================================
!    read cntrl vrbl in ensemble space (z)

      allocate(z(NENS_START:NENS))

      rewind icntrl
      READ(icntrl) z

      write(*,*) "INPUT: z min,max=",minval(z),maxval(z)

!------------------------------------------
        allocate(sum_loc(1:N_model))
        allocate(sumx(1:N_model))

                            call MATRIX_VECTOR &
   (Pa_file,namein,NENS_START,NENS,MPIRANK,NPROC,jsta,jend,N_model,z,sum_loc)

#ifdef MPI_USE
             CALL MPI_ALLREDUCE &
   (sum_loc,sumx,N_model,MPI_REAL4,MPI_SUM,MPI_COMM_WORLD,IERR)
#else
    sumx=sum_loc
#endif

       deallocate(sum_loc)
       deallocate(z)
       deallocate(jlen) ; deallocate(jdisp)

!========== write 1-d control variable ==============

!---- read background (first guess) vector

         if(MPIRANK.eq.0) then

      allocate(x_b(1:N_model))
      allocate(x(1:N_model))

      CLOSE(imodelb)
      OPEN(UNIT=imodelb,FILE='x1d_b',FORM='UNFORMATTED',IOSTAT=IER)
      IF(IER.NE.0) WRITE(*,*) ' imodelb OPEN UNIT ERROR IER=',IER
      read(imodelb) NN
      read(imodelb) x_b
      CLOSE(imodelb,status='KEEP')

      write(*,*) "BACKGROUND: x_b min,max=",minval(x_b),maxval(x_b)
      write(*,*) "DEL X: sumx min,max=",minval(sumx),maxval(sumx)

      x(:)=x_b(:) + sumx(:)

      write(*,*) "OUTPUT X: x min,max=",minval(x),maxval(x)
!-- write 1-d control variable
     CLOSE(imodel)
      OPEN(UNIT=imodel,FILE='x1d',FORM='UNFORMATTED',IOSTAT=IER)
      IF(IER.NE.0) WRITE(*,*) ' imodel OPEN UNIT ERROR IER=',IER
      write(imodel) NN
      write(imodel) x
      CLOSE(imodel,status='KEEP')

!---

      deallocate(x)
      deallocate(x_b)

         end if

      deallocate(sumx)

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

#ifdef MPI_USE
      CALL MPI_FINALIZE(IERR)
#endif

            stop
            end