program estimate_stddev
! **********************************************************************
! *              .      .    .
! * PROGRAM: estimate_stddev
! *
! * PRGMMR:     D. ZUPANSKI     ORG: CIRA/CSU   DATE: 2004-01-18
! *
! * ABSTRACT:   THIS SUBROUTINE ESTIMATES THE STANDARD DEVIATION
! *             (squre-root variance) FROM DIFFERENT model_history FILES
! *
! *
! * INPUT: 
! * model_history_yyyymmddhhmmss FILES
! *
! * OUTPUT:
! * std_dev_input FILE
! *
! * PROGRAM LOG:
! *
! * 08/14/2003    ..... M. ZUPANSKI, D. ZUPANSKI
! * 01/18/2004    ..... D. ZUPANSKI
! *
! -------------------------------------------------------------------------
! COMPILE:
!
! cp estimate_stddev.f  to a temporary directory  ($TMP)
! cd $TMP
! f90 -o test.x estimate_stddev.f
!
! RUN:
! cp ..../history/GEOS/model_true.cycle* .
! test.x
! cp  std_dev_input ..../ensda/namelists/.
! **********************************************************************

      integer,parameter :: ihist=25       ! input history file unit
      integer,parameter :: idev=105       ! output standard deviation
      real,dimension(:),allocatable::pt_mean,pt_var,pt_stdev
      real,dimension(:),allocatable::q_mean,q_var,q_stdev
      real,dimension(:,:),allocatable::pt_hist,q_hist

  
      character*20          :: filename
!------------------------------

      lm=55
      imax=10

      allocate(pt_hist(1:lm,1:imax))
      allocate(q_hist(1:lm,1:imax))

      allocate(pt_mean(1:lm))
      allocate(pt_var(1:lm))
      allocate(pt_stdev(1:lm))

      allocate(q_mean(1:lm))
      allocate(q_var(1:lm))
      allocate(q_stdev(1:lm))

      pt_mean(:)=0.
      pt_var(:)=0.
      pt_stdev(:)=0.
      q_mean(:)=0.
      q_var(:)=0.
      q_stdev(:)=0.

 100  format(E20.10)
               do i=1,imax
             if(i.lt.10) then
        write(filename,1000) i 
 1000 format('model_true.cycle',i1)
             else
        write(filename,1010) i
 1010 format('model_true.cycle',i2)
             endif
      write(*,*) " filename=",filename
      CLOSE(ihist)
      OPEN(UNIT=ihist,FILE=filename,FORM='FORMATTED',IOSTAT=IER)
      IF(IER.NE.0) WRITE(*,*) ihist,' OPEN UNIT ERROR IER=',IER
        read(ihist,100) PT_hist(1:lm,i)
        read(ihist,100) Q_hist(1:lm,i)
                 do l=1,lm
        pt_mean(l)= pt_mean(l)+PT_hist(l,i)
        q_mean(l)= q_mean(l)+q_hist(l,i)
        pt_var(l)=pt_var(l)+PT_hist(l,i)*PT_hist(l,i)
        q_var(l)=q_var(l)+Q_hist(l,i)*Q_hist(l,i)
                 enddo
               enddo

        pt_mean(:)=pt_mean(:)/float(imax)
        q_mean(:)=q_mean(:)/float(imax)
        pt_var(:)=pt_var(:)/float(imax-1)
        q_var(:)=q_var(:)/float(imax-1)

        write(*,*)" pt_mean=",pt_mean
        write(*,*)" pt_var=",pt_var
        write(*,*)" q_mean=",q_mean
        write(*,*)" q_var=",q_var

!                  do i=1,imax
!                   do l=1,lm
!       pt_stdev(l)=pt_stdev(l)+(PT_hist(l,i)-pt_mean(l))**2
!       q_stdev(l)=q_stdev(l)+(Q_hist(l,i)-q_mean(l))**2
!                   enddo
!                  enddo

!       pt_stdev(:)=sqrt(pt_stdev(:)/float(imax-1))
!       q_stdev(:)=sqrt(q_stdev(:)/float(imax-1))

        ff=float(imax)/float(imax-1)
        write(*,*)" ff=",ff
        pt_stdev(:)=sqrt(pt_var(:)-ff*pt_mean(:)**2)
        q_stdev(:)=sqrt(q_var(:)-ff*q_mean(:)**2)

        write(filename,1020) 
 1020 format('std_dev_input')
      write(*,*) " filename=",filename
      CLOSE(idev)
      OPEN(UNIT=idev,FILE=filename,FORM='FORMATTED',IOSTAT=IER)
      IF(IER.NE.0) WRITE(*,*) idev,' OPEN UNIT ERROR IER=',IER

        write(idev,*)"PT"
        write(idev,100) pt_stdev
        write(idev,*)"Q"
        write(idev,100) q_stdev

        write(*,*) " pt_stdev=",pt_stdev
        write(*,*) " q_stdev=",q_stdev

      write(*,*) " end estimate_stddev"

      end program estimate_stddev