PROGRAM PRIOR_STEP

! **********************************************************************
! *              .      .    .                                       
! * PROGRAM:	PRIOR_STEP
! * PRGMMR:	M. ZUPANSKI	ORG: CIRA/CSU	DATE: 2001-06-06
! *
! * ABSTRACT:	THIS PROGRAM PERFORMS THE PRIOR ERROR STEP CALCULATIONS
! *          	(FORECAST (BACKGROUND) ERROR AND MODEL ERROR)
! *
! * PROGRAM LOG:
! *
! * 06/06/2001    ..... M. ZUPANSKI: 
! * 02/10/2002    ..... M. ZUPANSKI: upgraded for ensembles 
! * 09/09/2003    ..... M. ZUPANSKI: Fortran-90 + further upgrade with latest code
! *
!-------------------------------------------------------------------
!   INPUT:    z(i)     = control variable (ensemble subspace)
!
!   INPUT:    sqrtIPA_inv = (I+A)**(-1/2) matrix
!
!   INPUT:    ALFA     = Old step-length, used for trial perturbation
!                       (ALFA is output in first, or restart minimization iteration)
!
!  OUTPUT:    SUM_UP, SUM_DOWN = sums needed for step-length update 
!
!-------------------------------------------------------------------
      integer,parameter::isha=20           ! restart and iteration parameters
      integer,parameter::IpA_file=21       ! (I+A)**(-1/2) file
      integer,parameter::icntrl=22         ! cntrl vrbl z (ensemble space)
      integer,parameter :: ialfa=23        ! step-length file
      integer,parameter :: idesc=24        ! descent direction file
      integer,parameter :: icost=51        ! step-length sums (inner-products)

      real ALFA
      real SGITOT
      integer ITERNO
      integer NENS
      integer NENS_START
      integer iout,ioutm1,index
      character*8 alfaname
      character*11 descname

      real,allocatable::sqrtIPA_inv(:,:)
      real,allocatable::z(:)
      real,allocatable::d(:)
      real,allocatable::s(:)
      real,allocatable::x(:)

      real,allocatable::gg(:)

!
!  DECLARE NAMELISTS
!
      NAMELIST /ENSEMBLE_SIZE/ NENS_START,NENS
      NAMELIST /CURRENT_ITERATION/ iout,ioutm1,index

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

!--- read restart parameters

      close(isha)
      OPEN(UNIT=isha,FILE='shanno',FORM='UNFORMATTED',IOSTAT=IER)
      IF(IER.NE.0)WRITE(*,*)' ISHA OPEN UNIT ERROR IER=',IER
      read(isha) SGITOT,ITERNO
      close(isha,status='KEEP')

!-- read current iteration

       REWIND 13
       READ(13,CURRENT_ITERATION)

       write(*,*) "prior_step: iout,ioutm1,index=",iout,ioutm1,index

!-- read ensemble size

       REWIND 15
       READ(15,ENSEMBLE_SIZE)

       write(*,*) "prior_step: NENS_START, NENS=",NENS_START,NENS

!--- read (I+A)**(-1/2) ---------------------------
      allocate(sqrtIPA_inv(NENS_START:NENS,NENS_START:NENS))

      REWIND IpA_file
      READ(IpA_file) ((sqrtIPA_inv(i,j),i=NENS_START,NENS),j=NENS_START,NENS)

!    read cntrl vrbl in ensemble space (z)

        allocate(z(NENS_START:NENS))

      CLOSE(icntrl)
      OPEN(UNIT=icntrl,FILE='z_cntrl',FORM='UNFORMATTED',IOSTAT=IER)
      IF(IER.NE.0) WRITE(*,*)' icntrl OPEN UNIT ERROR IER=',IER
      READ(icntrl) z
      CLOSE(icntrl,status='KEEP')

!--------------------------------------------------
 100  format('alfa_',i2.2)
 200  format('descent_',i2.2)

!-- read alfa from current iteration

      write(alfaname,100) iout
      close(ialfa)
      open(UNIT=ialfa,FILE=alfaname,FORM='UNFORMATTED',IOSTAT=IER)
      IF(IER.NE.0)WRITE(*,*)' ialfa OPEN UNIT ERROR IER=',IER
      read(ialfa) ALFA
      CLOSE(ialfa,status='KEEP')

!--- read descent direction

      write(descname,200) iout

      allocate(d(NENS_START:NENS))
      close(idesc)
      open(UNIT=idesc,FILE=descname,FORM='UNFORMATTED',IOSTAT=IER)
      IF(IER.NE.0) WRITE(*,*)' idesc OPEN UNIT ERROR IER=',IER
      read(idesc) (d(K),K=NENS_START,NENS)
      close(idesc,status='KEEP')

!!-- S = Z_trial - Z_old -------------
      allocate(s(NENS_START:NENS))
      s(:)=ALFA*d(:)

      deallocate(d)
!------------------------------------------
!!-- X = (I+A)**(-1)*(Z_trial - Z_old)
!-- gg=(I+A)**(-1/2)*z

         allocate(gg(NENS_START:NENS))

         do i=NENS_START,NENS
           gg(i)=0.0
          do j=NENS_START,NENS
           gg(i)=gg(i)+sqrtIPA_inv(i,j)*s(j)
          end do
         end do


!-- x=(I+A)**(-1/2)*gg=(I+A)**(-1)*z

         allocate(x(NENS_START:NENS))

         do i=NENS_START,NENS
           x(i)=0.0
          do j=NENS_START,NENS
           x(i)=x(i)+sqrtIPA_inv(i,j)*gg(j)
          end do
         end do

       deallocate(gg)

!=========================================================
!=====CALCULATE THE STEP_LENGTH CONTRIBUTIONS 
!========================================================

          BTOP=0.0
          BBOT=0.0

          do N=NENS_START,NENS
       BTOP=BTOP - z(N)*x(N)
       BBOT=BBOT + s(N)*x(N)
          end do

       deallocate(z)
!---- write step-length contributions

      close(icost)
      open(icost,file='backsum',status='unknown',form='unformatted')
      IF(IER.NE.0)WRITE(3,*)' icost OPEN UNIT ERROR IER=',IER
      write(icost) BTOP,BBOT
      close(icost,status='KEEP')

!------------------------------------------
       WRITE(*,1000) BTOP,BBOT
       WRITE(3,1000) BTOP,BBOT
 1000  FORMAT(' ','BTOP= ',E15.6,'  BBOT= ',E15.6/)
!------------------------------------------

      WRITE(3,*)"END PRIOR STEP CALCULATION "
      WRITE(*,*)"END PRIOR STEP CALCULATION "
!===========================================================
                      STOP
                      END