program grad_all_categ
!
!**********************************************************************
! *              .      .    .
! * PROGRAM:    grad_all
! * PRGMMR:     M. ZUPANSKI     ORG: CIRA/CSU   DATE: 2003-09-03
! *
! * ABSTRACT:   GRADIENT IN ENSEMBLE SPACE 
! *             Calculate total cost function gradient with respect to cntrl vrbl
! *             in ensemble subspace. 
! *             Combine background gradient with obs gradient
! *
! *
! * PROGRAM LOG:
! *
! * 09/03/2003    ..... M. ZUPANSKI:
! *
! **********************************************************************
!   INPUT:    z(i)     = control variable (ensemble subspace)
!
!   INPUT:    g_obs(i) = obs gradient (ensemble subspace)
!
!   INPUT:    sqrtIPA_inv = (I+A)**(-1/2) matrix
!                    
!  OUTPUT:    g(i) = (I+A)**(-1)*z - (I+A)**(-1/2)*g_obs
!
!--------------------------------------------------------------------
!-----
  integer,parameter::icntrl=21       ! cntrl vrbl (ensemble space)
  integer,parameter::igrad_obs=22    ! obs gradient
  integer,parameter::IpA_file=23     ! (I+A)**(-1/2)
  integer,parameter::gradout=51      ! total cost-function gradient  
!-----

 real,allocatable::sqrtIPA_inv(:,:)
 real,allocatable::z(:)
 real,allocatable::z_categ(:,:)
 real,allocatable::g(:)
 real,allocatable::g_categ(:,:)
 real,allocatable::gg(:)
 real,allocatable::g_obs(:)
 real,allocatable::gg_obs(:)
 real,allocatable::gg_obs_categ(:,:)
 real,allocatable::g_b(:)

 integer NENS,NALL
 integer NENS_START
 integer ncateg,nn

       CHARACTER*20 fileC

!
!      DECLARE NAMELIST
!
      NAMELIST /ENSEMBLE_SIZE/ NENS_START,NENS

!=========================================================
!==============start calculation===================
        write(*,*) "start grad_all"

!-- read ensemble size

       write(*,*) "grad_all: start reading input 15"

       REWIND 15
       READ(15,ENSEMBLE_SIZE)

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

      read(17,*) ncateg

      NALL=NENS-NENS_START+1
      N_start=NENS_START
      N_end=NALL*ncateg

 1140 FORMAT('sqrtIPA_inv_',I4.4)

        allocate(g(N_start:N_end))
        allocate(g_categ(NENS_START:NENS,1:ncateg))
        g(:)=0.0
        g_categ(:,:)=0.0

!    read the obs gradient (gg_obs)

        allocate(gg_obs(N_start:N_end))

        REWIND igrad_obs
        READ(igrad_obs) (gg_obs(k),k=N_start,N_end)

        allocate(gg_obs_categ(NENS_START:NENS,1:ncateg))

!    read the cntrl vrbl in ensemble space (z)

        allocate(z(N_start:N_end))

      CLOSE(icntrl)
      OPEN(UNIT=icntrl,FILE='zcntrl_file',FORM='UNFORMATTED',IOSTAT=IER)
      IF(IER.NE.0) WRITE(*,*)' icntrl OPEN UNIT ERROR IER=',IER
      READ(icntrl) z

        allocate(z_categ(NENS_START:NENS,1:ncateg))

      ij=N_start-1
             do j=1,ncateg
            do i=NENS_START,NENS
      ij=ij+1
      z_categ(i,j)=z(ij)
      gg_obs_categ(i,j)=gg_obs(ij)
            end do
             end do

        deallocate(gg_obs)
        deallocate(z)

               do nn=1,ncateg
!--- read (I+A)**(-1/2) ---------------------------
      allocate(sqrtIPA_inv(NENS_START:NENS,NENS_START:NENS))

      WRITE(fileC,1140) nn
      CLOSE(IpA_file)
      OPEN(IpA_file,file=fileC,status='unknown',form='unformatted')
      READ(IpA_file) ((sqrtIPA_inv(i,j),i=NENS_START,NENS),j=NENS_START,NENS)

!==================================================
!==== prior (background) gradient =================
! (note that (I+A)**(-1/2) is symmetric => does not matter if transpose or not)

         allocate(gg(NENS_START:NENS))

!-- gg=(I+A)**(-1/2)*z
         do i=NENS_START,NENS
           gg(i)=0.0
          do j=NENS_START,NENS
           gg(i)=gg(i)+sqrtIPA_inv(i,j)*z_categ(j,nn)
          end do
         end do

         allocate(g_b(NENS_START:NENS))

!-- g_b=(I+A)**(-1/2)*gg
         do i=NENS_START,NENS
           g_b(i)=0.0
          do j=NENS_START,NENS
           g_b(i)=g_b(i)+sqrtIPA_inv(i,j)*gg(j)
          end do
         end do

       deallocate(gg)

!=========================================================
!  obs gradient
!-- g_obs=(I+A)**(-1/2)*gg_obs

         allocate(g_obs(NENS_START:NENS))

!        do i=NENS_START,NENS
!          g_obs(i)=0.0
!         do j=NENS_START,NENS
!          g_obs(i)=g_obs(i)+sqrtIPA_inv(i,j)*gg_obs(j)
!         end do
!        end do

      g_obs(:)=gg_obs_categ(:,nn)

     deallocate(sqrtIPA_inv)
!====================================================================
!  TOTAL GRADIENT g IN ENSEMBLE SPACE
!====================================================================

          do i=NENS_START,NENS
           g_categ(i,nn)=g_b(i)+g_obs(i)
          end do

       deallocate(g_b) ; deallocate(g_obs)

               end do  !!!do nn=1,ncateg

!====================================================================

          do jj=NENS_START,min(20,NENS)
          ii=min(5,ncateg)
          write(*,101) jj,g_categ(jj,ii)
          end do

! write out the total gradient (g)

      ij=N_start-1
             do j=1,ncateg
            do i=NENS_START,NENS
      ij=ij+1
      g(ij)=g_categ(i,j)
            end do
             end do

        rewind gradout
        write(gradout) (g(k),k=N_start,N_end)

       deallocate(z_categ) ; deallocate(g) ; deallocate(g_categ)


 101    format('g(',I4,')=',E17.5)
!====================================================================

 stop
 end