PROGRAM INFLATE

! **********************************************************************
! *              .      .    .                                       
! * PROGRAM:	ERROR COVARIANCE 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:    OLD_FACTOR     = previous inflation factor
!
!   INPUT:    res            = normalized observation residual 
!
!   INPUT:    xlambda        = eigenvalues of matrix A
!
!  OUTPUT:    FACTOR         = updated inflation factor
!
!-------------------------------------------------------------------
      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
      integer,parameter :: isqrtPx=52     ! square-root analysis error covariance

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

      real P_FACTOR,OLD_FACTOR,FACTOR
      integer NENS
      integer NENS_START
      integer mobs
      integer NALL
      integer icycle
      integer model_dim
      character*11  Pxname

      real,dimension(:),allocatable :: xlambda
      real,dimension(:),allocatable :: res
      real,dimension(:),allocatable :: sqrtPx_in,sqrtPx

!
!  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(ires)
      OPEN(ires,file='res_cntrl',status='unknown',form='unformatted')
      read(ires) mobs
      allocate(res(max(1,mobs)))
      read(ires) (res(k),k=1,mobs)
      CLOSE(ires,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='inflation_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='inflation_out',status='unknown',form='unformatted')
      write(infla) P_FACTOR
      CLOSE(infla,STATUS='KEEP')
!--------------------------------------------------
 2000 format('sqrtPa_',i4.4)

      do N=NENS_START,NENS

      WRITE(Pxname,2000)  N
      CLOSE(isqrtPx)
      OPEN(isqrtPx,file=Pxname,status='unknown',form='unformatted')
      READ(isqrtPx) model_dim
      allocate(sqrtPx(1:model_dim))
      READ(isqrtPx) sqrtPx
      CLOSE(isqrtPx,status='KEEP')

      write(*,*) "BEFORE: N=",N," sqrtPa min,max=",minval(sqrtPx),maxval(sqrtPx)

           do i=1,model_dim
           sqrtPx(i) = P_FACTOR * sqrtPx(i)
           end do

      write(*,*) "AFTER: N=",N," sqrtPa min,max=",minval(sqrtPx),maxval(sqrtPx)

      CLOSE(isqrtPx)
      OPEN(isqrtPx,file=Pxname,status='unknown',form='unformatted')
      WRITE(isqrtPx) model_dim
      WRITE(isqrtPx) sqrtPx
      CLOSE(isqrtPx,status='KEEP')

      deallocate(sqrtPx)

      end do
!--------------------------------------------------


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