program gauss_histogram
!**********************************************************************
! *              .      .    .
! * PROGRAM:    GAUSS_HISTOGRAM
! * PRGMMR:     M. ZUPANSKI     ORG: CIRA/CSU   DATE: 2003-09-09
! *
! * ABSTRACT:   CALCULATE GAUSSIAN PDF HISTOGRAM
! *
! * PROGRAM LOG:
! *
! * 09/09/2003    ..... M. ZUPANSKI: Original 'transform.F'
! * 10/17/2003    ..... M. ZUPANSKI: Innovation normalization
! * 10/22/2003    ..... M. ZUPANSKI: Innovation histogram
! * 10/27/2003    ..... M. ZUPANSKI: Gaussian PDF histogram
! *
! **********************************************************************
!----------------------------------
!  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::hist_file=51      ! Innovation histogram  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

     character*30 filename

     real,dimension(:),allocatable::bin
     real,dimension(:),allocatable::sbin
     real,dimension(:),allocatable::sigma
     real,dimension(:),allocatable::x_innov
     integer,dimension(:),allocatable::ibin
!-----
!
!      DECLARE NAMELIST
!
     NAMELIST /INNOV_BINS/ bin_increment,bin_range

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

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

      i1max=bin_range/bin_increment
      allocate(ibin(-i1max:i1max))
      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

      
      allocate(sigma(1:3))
      sigma(1)=1.0
      sigma(2)=sqrt(2.0)
      sigma(3)=0.5 

       do ii=1,3

      const_pi=1./(sigma(ii)*sqrt(2.*3.14))
      sigma2=sigma(ii)**2
      write(*,*) ii," sigma=",sigma(ii)," sigma**2=",sigma2
      write(*,*) ii," const_pi=",const_pi

      allocate(sbin(-i1max:i1max))
      do i=-i1max,i1max
      sbin(i)=const_pi*exp(-bin(i)**2/(2.*sigma2))
      end do

      sumbin=0.0
      do i=-i1max,i1max
      write(*,*) i," sbin=",sbin(i)
      sumbin=sumbin+sbin(i)
      end do
      write(*,*) ii," sumbin=",sumbin
        
      write(filename,110) ii
      CLOSE(hist_file)
      OPEN(UNIT=hist_file,FILE=filename,FORM='FORMATTED',IOSTAT=IER)
      IF(IER.NE.0) WRITE(*,*) hist_file,' OPEN UNIT ERROR IER=',IER
      do i=-i1max,i1max
      write(hist_file,200) sbin(i)
      end do
      CLOSE(hist_file,status='KEEP')

      deallocate(sbin)

       end do   !! ii-loop

 200  format(E20.10)
 110  format('bin_hist_gauss_',i1.1)


      write(*,*) " sumbin=",sumbin," ntotal=",ntotal

            end program gauss_histogram