program ibin_cycle_sum
!**********************************************************************
! *              .      .    .
! * PROGRAM:    IBIN_CYCLE_SUM
! * PRGMMR:     M. ZUPANSKI     ORG: CIRA/CSU   DATE: 2003-09-09
! *
! * ABSTRACT:   CALCULATE IBINS FOR CURRENT CYCLE
! *
! * PROGRAM LOG:
! *
! * 09/09/2003    ..... M. ZUPANSKI: Original 'transform.F'
! * 10/17/2003    ..... M. ZUPANSKI: Innovation normalization
! * 10/22/2003    ..... M. ZUPANSKI: Innovation histogram
! * 11/10/2004    ..... M. ZUPANSKI: Ibins by cycles
! *
! **********************************************************************
!----------------------------------
!  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
!
!----------------------------------
!-----
     integer,parameter::innov_file=21     ! Normalized innovation  file #
     integer,parameter::ibin_file=31      ! Updated ibin file
     integer,parameter::ix_sum=41         ! Updated x_sum file
!-----
     real    :: bin_increment,bin_range
     real    :: bin_left,bin_right
     real    :: bin_left_max,bin_right_max
     integer :: N_obs
     integer :: i1max,i,N

     real    :: x1_sum,x2_sum
     integer :: Nx

     real,dimension(:),allocatable::bin
     real,dimension(:),allocatable::x_innov
     integer,dimension(:),allocatable::ibin
     integer,dimension(:),allocatable::ibin_old
     character*1 categ
!-----
!
!      DECLARE NAMELIST
!
     NAMELIST /INNOV_BINS/ bin_increment,bin_range

!==============start calculation===================
      write(*,*) "START IBIN_CYCLE_SUM"

!---- define bins (from dbin, binmax)
       REWIND 16
       READ(16,INNOV_BINS)

      i1max=bin_range/bin_increment
      allocate(bin(-i1max:i1max))

      bin(0)=0.0
      do i=1,i1max
      bin(i)=bin(i-1)+bin_increment
      bin(-i)=bin(-i+1)-bin_increment
      end do

      write(*,*) "bin_increment,bin_range,i1max=",bin_increment,bin_range,i1max

      do i=-i1max,i1max
      write(*,*) i," bin=",bin(i)
      end do

 100  format('innov.',i3.3)

      bin_left_max=bin(-i1max)-bin_increment/2.
      bin_right_max=bin(i1max)+bin_increment/2.

      CLOSE(innov_file)
      OPEN(UNIT=innov_file,FILE='innov',FORM='UNFORMATTED',IOSTAT=IER)
      IF(IER.NE.0) WRITE(*,*) innov_file,' OPEN UNIT ERROR IER=',IER
      read(innov_file) N_obs
      allocate(x_innov(1:N_obs))
      read(innov_file) x_innov
      CLOSE(innov_file,status='KEEP')

      write(*,*) N," Innov min,max=",minval(x_innov),maxval(x_innov)

      allocate(ibin_old(-i1max:i1max))
      allocate(ibin(-i1max:i1max))
      ibin(:)=0

!       call getenv("icateg",categ)
!       if(categ.eq.'1') then
           ibin_old(:)=0
           x1_sum=0.0
           x2_sum=0.0
           Nx=0
!       else
           CLOSE(ibin_file)
           OPEN(UNIT=ibin_file,FILE='ibin_sum_file',FORM='FORMATTED',IOSTAT=IER)
           IF(IER.NE.0) WRITE(*,*) ibin_file,' OPEN UNIT ERROR IER=',IER
           read(ibin_file,150,end=220) (ibin_old(i),i=-i1max,i1max)
           CLOSE(ibin_file,status='KEEP')

           CLOSE(ix_sum)
           OPEN(UNIT=ix_sum,FILE='x_sum_file',FORM='FORMATTED',IOSTAT=IER)
           IF(IER.NE.0) WRITE(*,*) ix_sum,' OPEN UNIT ERROR IER=',IER
           read(ix_sum,*,end=220) Nx,x1_sum,x2_sum
           CLOSE(ix_sum,status='KEEP')
!       end if

 220    continue

      do i=1,N_obs

      Nx=Nx+1
      x1_sum=x1_sum+x_innov(i)
      x2_sum=x2_sum+x_innov(i)**2

        do k=-i1max,i1max
      bin_left=bin(k)-bin_increment/2.
      bin_right=bin(k)+bin_increment/2.
       if(x_innov(i).ge.bin_left.and.x_innov(i).lt.bin_right) then
       ibin(k)=ibin(k)+1
       end if
        end do  !! bin-loop
        if(x_innov(i).lt.bin_left_max) then
        ibin(-i1max)=ibin(-i1max)+1
          elseif(x_innov(i).ge.bin_right_max) then
          ibin(i1max)=ibin(i1max)+1
        end if

       end do    !! obs-loop
        
      deallocate(x_innov)

      do i=-i1max,i1max
      ibin(i)=ibin(i)+ibin_old(i)
      write(*,*) i," ibin_old=",ibin_old(i)," ibin=",ibin(i)
      end do

      CLOSE(ibin_file)
      OPEN(UNIT=ibin_file,FILE='ibin_sum_file',FORM='FORMATTED',IOSTAT=IER)
      IF(IER.NE.0) WRITE(*,*) ibin_file,' OPEN UNIT ERROR IER=',IER
      write(ibin_file,150) (ibin(i),i=-i1max,i1max)
      CLOSE(ibin_file,status='KEEP')
 150  format(I10)

           CLOSE(ix_sum)
           OPEN(UNIT=ix_sum,FILE='x_sum_file',FORM='FORMATTED',IOSTAT=IER)
           IF(IER.NE.0) WRITE(*,*) ix_sum,' OPEN UNIT ERROR IER=',IER
           write(ix_sum,*) Nx,x1_sum,x2_sum
           CLOSE(ix_sum,status='KEEP')
        
      write(*,*) "END IBIN_CYCLE_SUM"

            end program ibin_cycle_sum