program analysis_cov
!**********************************************************************
! *              .      .    .
! * PROGRAM:    ANALYSIS_COV
! * PRGMMR:     M. ZUPANSKI     ORG: CIRA/CSU   DATE: 2003-08-28
! *
! * ABSTRACT:   SQUARE-ROOT ANALYSIS ERROR COVARIANCE CALCULATION
! *
! * PROGRAM LOG:
! *
! * 08/28/2003    ..... M. ZUPANSKI:
! *
! **********************************************************************
!----------------------------------
!  Hessian preconditioning matrix
!
!  INPUT: 
!         sqrt(Pf) - fcst_xxxx
!         sqrt(I+A)**(-1/2) - sqrtIPA_inv
!
!  OUTPUT:
!         sqrtPa=sqrtPf*sqrtIPA_inv
!
!         (Pa)**(1/2)=(Pf)**(1/2) * (I+A)**(-T/2)
!
!----------------------------------

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

     integer,parameter::IpA_file=21       ! (I+A)**(-T/2) file #
     integer,parameter::Pf_file=22        ! ensemble fcst ptrb file #
     integer,parameter::Pa_file=51        ! sqrtPa file #
!-----
     integer :: NENS
     integer :: NENS_START
     integer :: NALL
     integer :: N_LOC,jsta,jend
     integer :: NPROC,MPIRANK
     character*20 :: name_in,name_out
     character*20 :: name_ens
     real,allocatable :: sqrtIPA_inv(:,:)
     integer,allocatable :: jlen(:),jdisp(:)
!-----
!
!      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
!================================================

             if(MPIRANK.eq.0) then
       write(*,*) "analysis_cov: start reading ensemble size (file 15)"
             end if
       REWIND 15
       READ(15,ENSEMBLE_SIZE)

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

!-- Total number of degrees of freedom

        NALL=NENS-NENS_START+1

!====================================================================
!===  create a column matrix Pa**(1/2) ==
!====================================================================

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

!-- input filename
      write(name_in,101)
 101  format('sqrtPf_')
!-- output filename
      write(name_out,102)
 102  format('sqrtPa_')
!-- ensemble matrix name
      write(name_ens,103)
 103  format('sqrtIpA_inv')

!--- read (I+A)**(-1/2) ---------------------------
      allocate(sqrtIPA_inv(NENS_START:NENS,NENS_START:NENS))
         if(MPIRANK.eq.0) then
         write(*,*) "read sqrt(I+A)**(-1/2) matrix from file ",IpA_file
         end if
      REWIND IpA_file
      READ(IpA_file) ((sqrtIPA_inv(i,j),i=NENS_START,NENS),j=NENS_START,NENS)


                            call MATRIX_MATRIX &
   (Pf_file,Pa_file,name_in,name_out,NENS_START,NENS,MPIRANK,NPROC,jsta,N_LOC,sqrtIPA_inv)


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

      deallocate(jlen) ; deallocate(jdisp)
      deallocate(sqrtIPA_inv)

#ifdef MPI_USE
      CALL MPI_FINALIZE(IERR)
#endif

            stop
            end