PROGRAM INFLATE

! **********************************************************************
! *              .      .    .                                       
! * PROGRAM:	INFLATION
! * PRGMMR:	M. ZUPANSKI	ORG: CIRA/CSU	DATE: 2004-04-14
! *
! * ABSTRACT:	INFLATION FACTOR FOR USE IN ETKF (Wang and Bishop, 2003, JAS)
! *
! * PROGRAM LOG:
! *
! * 04/15/2004    ..... M. ZUPANSKI: 
! *
!-------------------------------------------------------------------
!   INPUT:    z(i)     = control variable (ensemble subspace)
!
!   INPUT:    ALFA     = Step-length
!
!   INPUT:    d(i)     = descent direction
!
!  OUTPUT:    z_step   = z + ALFA * d
!
!-------------------------------------------------------------------
      integer,parameter:: infla_old=21    ! previous inflation factor
      integer,parameter :: ires=22        ! normalized observation residuals
      integer,parameter :: ieigen=23      ! eigenvalues of matrix A
      integer,parameter :: infla=51       ! updated inflation factor

!---------------------------------------------------------------------------------------------

      real P_FACTOR,OLD_FACTOR,FACTOR
      integer NENS
      integer NENS_START
      integer mobs
      integer NALL
      integer icycle

      real,allocatable::xlambda(:)
      real,allocatable::res(:)

!
!  DECLARE NAMELISTS
!
      NAMELIST /ENSEMBLE_SIZE/ NENS_START,NENS
      NAMELIST /CURRENT_CYCLE/ icycle

!==============start calculation===================
      WRITE(3,*)"START COVARIANCE INFLATION FACTOR "
      WRITE(*,*)"START COVARIANCE INFLATION FACTOR "

!-- read current cycle
      write(*,*) "read current cycle"
      iparm=152
      CLOSE(iparm)
      OPEN(UNIT=iparm,FILE='cycle.name',FORM='FORMATTED',IOSTAT=IER)
      IF(IER.NE.0)WRITE(*,*) iparm,' OPEN UNIT ERROR IER=',IER
      READ(iparm,CURRENT_CYCLE)

       write(*,*) "INFLATION: current cycle=",icycle

!-- read ensemble size

       REWIND 15
       READ(15,ENSEMBLE_SIZE)

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

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

!
! read in the obs incr 
!

      CLOSE(idim)
      OPEN(idim,file='res',status='unknown',form='unformatted')
      read(idim) mobs
      allocate(res(max(1,mobs)))
      read(idim) (res(k),k=1,mobs)
      CLOSE(idim,STATUS='KEEP')

!-- A_eigenfile --------------------
!   read eigenvalues of A
!! (NOTE: (I+Lambda)**(-1/2) is on th file. Need to convert to Lambda !!!!)

      allocate(xlambda(NENS_START:NENS))

      CLOSE(ieigen)
      OPEN(UNIT=ieigen,FILE='A_eigenfile',FORM='UNFORMATTED',IOSTAT=IER)
      IF(IER.NE.0) WRITE(*,*) ieigen,' OPEN UNIT ERROR IER=',IER
      READ(ieigen) NN
         if(NN.ne.NALL) then
         write(*,*) "Innov_norml DIMENSION PROBLEM: NN=",NN," NALL=",NALL
         stop
         end if
      READ(ieigen) (xlambda(i),i=NENS_START,NENS)
      READ(ieigen) 
      CLOSE(ieigen,STATUS='KEEP')

      do N=NENS_START,NENS
       diff=1./xlambda(N)**2-1.
       xlambda(N)=max(diff,0.0)
      end do

!-- 
        sum_lambda=0.0
        do i=NENS_START,NENS-1
        sum_lambda=sum_lambda+xlambda(i)
        end do
     deallocate(xlambda)
!-- 
        sum_res=0.0
        do i=1,mobs
        sum_res=sum_res+res(i)*res(i)
        end do
     deallocate(res)

        write(*,*) " sum_lambda=",sum_lambda
        write(*,*) " sum_res=",sum_res," N obs=",mobs

         if(sum_lambda.gt.0.0) then
        FACTOR=(sum_res - float(mobs))/sum_lambda
        write(*,*) "FACTOR=",FACTOR
         else
        write(*,*) "PROBLEM with sum_lambda=",sum_lambda
        stop
         end if

!
! Inflation factor update
!
       if(icycle.gt.1) then
      CLOSE(infla_old)
      OPEN(infla_old,file='p_factor',status='unknown',form='unformatted')
      read(infla_old) OLD_FACTOR
      CLOSE(infla_old,STATUS='KEEP')
       else
      OLD_FACTOR=1.0
       end if

      write(*,*) "OLD_FACTOR=",OLD_FACTOR

       if(FACTOR.ge.0.0) then
      P_FACTOR=OLD_FACTOR*sqrt(FACTOR)
        write(*,*) "UPDATED P_FACTOR=",P_FACTOR
       else
      write(*,*) "PROBLEM with FACTOR=",FACTOR
      P_FACTOR=OLD_FACTOR*sqrt(abs(FACTOR))
       end if

! write updated inflation factor
      CLOSE(infla)
      OPEN(infla,file='p_out',status='unknown',form='unformatted')
      write(infla) P_FACTOR
      CLOSE(infla,STATUS='KEEP')
!--------------------------------------------------


      WRITE(3,*)"END COVARIANCE INFLATION FACTOR "
      WRITE(*,*)"END COVARIANCE INFLATION FACTOR "
!===========================================================
                      STOP
                      END