!
      module m_obs_pert
!
! Module to perturb observations with random perturbations
! For use with observation simulation software for the NASA version of the 
! ECMWF/NCEP/NASA/etc. OSSE system.
!
! Initial code provided by Ronald Errico (Jan. 2008) 
!
      use m_kinds
      use m_relhum
!
      private
      public pert_find_itype
      public pert_obs
!
      contains
!
!
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
!
      subroutine pert_eigen (nk,matrix,evects,evalues)
!
! Call library routine to compute eigenvalues and vectors of a symmetric 
! matrix (stored as an array of 8-byte values).
!
      integer,  intent(in)   :: nk
      real(rkind3), intent(in)   :: matrix(nk,nk)
      real(rkind3), intent(out)  :: evects(nk,nk)
      real(rkind3), intent(out)  :: evalues(nk)
!
      integer(4)     :: nwork1  
      integer(4)     :: info 
      integer(4)     :: nk4
      real(rkind3)    :: work1(nk*nk*2)
!
      nk4=nk
      nwork1=nk4*nk4*2
      evects=matrix
      call dsyev ( 'v','u',nk4,evects,nk4,evalues,work1,nwork1,info)
      if (info /= 0) then
        print *,'   '
        print *,' * * * * * * * * * * * * * * * * * * * * * * *  '
        print *,' info from dsyev =  ',info
        print *,'   '
      endif
!
      end subroutine pert_eigen
!
!
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
!
      Subroutine pert_find_itype (n_err3,itype,err_itype,err_index)
      implicit none
      integer :: n_err3
      integer :: itype
      integer :: err_itype(n_err3)
      integer :: err_index
      integer :: k  
!     
      err_index=0
      do k=1,n_err3
        if (itype == err_itype(k)) then
          err_index=k
        endif
      enddo
!
      end subroutine pert_find_itype 
!
!
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
!
      Subroutine pert_fix_evalues (nk,evalues)
!
!  Set small (relative to largest values) or negative values to 
!  a small value positive value
!
      integer,  intent(in)    :: nk
      real(rkind3), intent(inout) :: evalues(nk)
!
      integer  :: k
      real(rkind3) :: emax, emin
!
      emax=0.
      do k=1,nk
        if (abs(evalues(k)) > emax) emax=abs(evalues(k))
      enddo
!
      if (rkind3==8) then
        emin=emax*1.e-12
      else
        emin=emax*1.e-6
      endif
!
      do k=1,nk
        if (abs(evalues(k)) < emin) evalues(k)=emin
      enddo
!
      end subroutine pert_fix_evalues
!
!
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
!
      subroutine pert_gauss (mean,stdv,draw1)
!
! Retrieve a random number for a Gaussian distribution with input 
! mean and standard deviation.  Algorithm uses central limit theorem.
!
      integer, parameter :: ndraw=12
      real(rkind3), intent(in)   :: mean
      real(rkind3), intent(in)   :: stdv
      real(rkind3), intent(out)  :: draw1
!
      integer      :: n
      real(rkind3) :: x, xsum
!
      xsum=zero_k3
      do n=1,ndraw
        call random_number (x)
        xsum=xsum+x
      enddo
! The following x is normal(0,1)       
      x=(xsum-0.5_rkind3*ndraw)*sqrt(12._rkind3/ndraw)
      draw1=x*stdv+mean
!
      end subroutine pert_gauss
!
!
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
!
      subroutine pert_interp_stdv (nptab,p_tab,stdv_tab,p,stdv) 
!  
! Interpolate stdv from table defined on p-levels defined by p_tab.
!
      integer :: nptab
      real(rkind3) :: p_tab(nptab)
      real(rkind3) :: stdv_tab(nptab)
      real(rkind3) :: p
      real(rkind3) :: stdv
      real(rkind3) :: r1,r2
      integer :: k, klev
!
! first find klev as p-level just below p (i.e., next largers p_lev)

      if (p > p_tab(1) ) then 
        stdv=stdv_tab(1)
      elseif (p < p_tab(nptab) ) then
        stdv=stdv_tab(nptab)
      else
        klev=1
        do k=2,nptab-1
          if (p > p_tab(k)) exit
          klev=k
        enddo 
!
        r2=(p_tab(klev)-p)/(p_tab(klev)-p_tab(klev+1))
        r1=one_k1-r2
        stdv=r1*stdv_tab(klev)+r2*stdv_tab(klev+1)
!       
      endif
!
      end subroutine pert_interp_stdv 
!
!
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
!
      Subroutine pert_normalize_evects (nk,evects)
!
!  Normalize each eigenvector such that the sum of its squared 
!  components is 1
!
      integer,  intent(in)    :: nk
      real(rkind3), intent(inout) :: evects(nk,nk)
!
      integer  :: k, i
      real(rkind3) :: sum  
     
      do k=1,nk
        sum=zero_k3
        do i=1,nk
          sum=sum+evects(i,k)*evects(i,k)
        enddo
        sum=1./sqrt(sum)
        evects(:,k)=evects(:,k)*sum
      enddo
!
      end subroutine pert_normalize_evects
!
!
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
!
       Subroutine pert_obs (n_err1,n_err2,obs_ndim1,obs_ndim2,nobs,   &
                      obs_nlevs,obs_nfields,icolumns,bmiss,pert_fac,  &
                      err_tab,obs_levels,obs_values,corr_dist,        &
                      ltest,lrelhum,lrad)
!
!  Add perturbation to observations.
!  Perturbations are drawn from a Gaussian distribution
!  For observations at multiple levels (e.g., for RAOBS), the perts are 
!  drawn independentaly for each eigenvector of the covariance matrix, 
!  assuming a vertical correlation function that has shape that is 
!  approximately Gaussian in vertical separation.  No correlation is assumed 
!  between different fields of the same observation (i.e., between t,q or
!  u,v.
!
!  Specifically, the verical correlation function is c=e(a*z**2) where 
!  z= log(p1/p2) for obs at p-levels p1 and p2, 
!  and a=log(xfac_r)/log(xfac_p)**2;
!  such that when either p1/p2 or p2/p1 = xfac_p, then c=xfac_r
!
!  Standard deviations of errors of conventional observations are determined
!  by interpolated from tabular values to the p-level of the observation.
!  For radiances, no interpolation between channels is required to 
!  determine the standard deviations of the corresponding observation 
!  values.  Therefore the standard deviations are simply copied from the
!  table provided.
!
      implicit none
!
      integer :: n_err1
      integer :: n_err2
      integer :: obs_ndim1 
      integer :: obs_ndim2 
      integer :: obs_nlevs
      integer :: obs_nfields
      integer :: icolumns(obs_nfields)
      integer :: nobs
      real(rkind3) :: pert_fac
      real(rkind3) :: bmiss
      real(rkind3) :: err_tab(n_err1,n_err2)
      real(rkind3) :: obs_levels(obs_ndim1)
      real(rkind3) :: obs_values(obs_ndim1,obs_ndim2)
      real(rkind3) :: corr_dist(obs_nfields)
      logical :: ltest
      logical :: lrelhum
      logical :: lrad
!
      integer      :: i, j, nf, icol
      real(rkind3) :: corr_d
      real(rkind3) :: cov(obs_nlevs,obs_nlevs)
      real(rkind3) :: evalues(obs_nlevs)
      real(rkind3) :: evects(obs_nlevs,obs_nlevs)
      real(rkind3) :: obs_old(obs_nlevs)
      real(rkind3) :: pert
      real(rkind3) :: sqrtev
      real(rkind3) :: stdv(obs_nlevs)
      real(rkind3) :: x
      real(rkind3) :: xfac
      real(rkind3), parameter :: xfac_r=0.1_rkind3
!
! Loop over obs fields
      do nf=1,obs_nfields
!
! save original values for later comparison and to save q for bmiss
        obs_old(1:obs_nlevs)=obs_values(1:obs_nlevs,nf)
!    
! If obs. field #2 is for q, perturb as relative humidity
! Therefore convert specific humidity obs to relative humidity obs
        if (lrelhum .and. (nf==2)) then
          call pert_q2rh (obs_nlevs,obs_values(1:obs_nlevs,1),     &
                obs_levels(1:obs_nlevs),obs_values(1:obs_nlevs,2) )
        endif 
!
! Fill arrays of standard deviations for each obs error 
        icol=icolumns(nf)
        do i=1,obs_nlevs
          if (lrad) then   ! copy values for corresponding channels
            stdv(i)=err_tab(i,icol) 
          else             ! interpolate between p-levels in table
            call pert_interp_stdv (n_err1,err_tab(:,1), &
                  err_tab(:,icol),obs_levels(i),stdv(i))
          endif
        enddo
!
! Determine if obs are not vertically correlated
        if ( (obs_nlevs == 1) .or. (corr_dist(nf) <= one_k3) ) then 
!          
          do i=1,obs_nlevs 
            if (obs_values(i,nf) >= bmiss) then
              obs_values(i,nf)=bmiss
            else
              call pert_gauss (zero_k3,stdv(i),pert)
              pert=pert*pert_fac
              obs_values(i,nf)=obs_values(i,nf)+pert
            endif
          enddo
!
        else  ! obs errors are correlated in the vertical 
!
! Fill symmetric covariance matrix
          xfac=log(xfac_r)/(log(corr_dist(nf)))**2 
          do i=1,obs_nlevs
            do j=i,obs_nlevs
              x=log(obs_levels(i)/obs_levels(j))
              cov(i,j)=stdv(i)*stdv(j)*exp(xfac*x*x)
              cov(j,i)=cov(i,j)
            enddo
          enddo
!
! Compute eigen values and normalized eigen vectors of cov matrix
! Also, do not allow relatively small or negative eigenvalues
          call pert_eigen (obs_nlevs,cov,evects,evalues)
          call pert_normalize_evects (obs_nlevs,evects)
          call pert_fix_evalues (obs_nlevs,evalues)
!
! Perturb each eigenvector independently and sum result
          do i=1,obs_nlevs
            sqrtev=sqrt(evalues(i))
            call pert_gauss (zero_k3,sqrtev,pert)
            pert=pert*pert_fac
            do j=1,obs_nlevs
              obs_values(j,nf)=obs_values(j,nf)+pert*evects(j,i)
            enddo
          enddo
!
        endif      ! test on whether vertically correlated    
!
! If obs. field #2 is for q, perturb as relative humidity
! Therefore convert back to specific humdity here
        if (lrelhum .and. (nf==2)) then
          call pert_rh2q (obs_nlevs,obs_values(1:obs_nlevs,1), &
                        obs_levels(1:obs_nlevs),obs_values(1:obs_nlevs,2) )
        endif 
!
! Restore missing values in original data
        do i=1,obs_nlevs
          if (obs_old(i) >= bmiss) then
            obs_values(i,nf)=bmiss
          endif
        enddo  
! 
! Print some arrays if testing
        if (ltest) then
          print *,' '
          print *,' obs number=',nobs,' nf=',nf,' corr_dist=',corr_dist(nf)
          print *,' Observations: pressure, old value,', &
                   ' new value, diff, stdv'
          do i=1,obs_nlevs
            x=obs_values(i,nf)-obs_old(i)
            print ('(i2,1p5e15.4)'),i,obs_levels(i),       &
                     obs_old(i),obs_values(i,nf),x,stdv(i)
          enddo
        endif    ! test on ltest
 
      enddo      ! loop over nfields
!
      end subroutine pert_obs 
!
!
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
!
      subroutine pert_print_matrix (n1,n2,np,matrix,cname)
!
! Print matrix
!
      integer,  intent(in)   :: n1, n2, np
      real(rkind3), intent(in)    :: matrix(n1,n2)
      character(len=*)       :: cname
!
      integer     :: n  
!
      print *,cname,' np,n2=',np,n2
      do n=1,n2
        print ('(1p10e12.2)'),matrix(1:np,n)
      enddo
!
      end subroutine pert_print_matrix
!
!
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
!
      subroutine pert_print_vector (n1,vector,cname)
!
! Print vector
!
      integer,  intent(in)   :: n1
      real(rkind3), intent(in)    :: vector(n1)
      character(len=*)       :: cname
!
      print *,cname,' n1=',n1
      print ('(1p10e12.2)'),vector(:)
!
      end subroutine pert_print_vector
!
!
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
!
      subroutine pert_q2rh (nlevs,t,p,q)
!
! Convert q to rh
!
      integer,      intent(in)    :: nlevs
      real(rkind3), intent(in)    :: t(nlevs), p(nlevs)
      real(rkind3), intent(inout) :: q(nlevs)
!
      integer      :: k, ier
      real(rkind2) :: t2, p2, q2
!    
      do k=1,nlevs  
        ier=0
        t2=t(k)      
        p2=p(k)      
        q2=q(k)
        call relhum_q2rh (t2,p2,q2,ier)
        if (ier == 0 ) then
          q(k)=q2
        else   
          q(k)=1.e-4
        endif
      enddo 
!
      end subroutine pert_q2rh
!
!
!  X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X
!
      subroutine pert_rh2q (nlevs,t,p,q)
!
! Convert rh to q and ensure positive value
!
      integer,      intent(in)    :: nlevs
      real(rkind3), intent(in)    :: t(nlevs), p(nlevs)
      real(rkind3), intent(inout) :: q(nlevs)
!
      integer      :: k, ier
      real(rkind1) :: t2, p2, q2
!    
      do k=1,nlevs  
        ier=0
        t2=t(k)      
        p2=p(k)      
        q2=q(k)
        call relhum_rh2q (t2,p2,q2,ier)
        if ( (ier == 0 ) .and. (q2 > zero_k3) ) then
          q(k)=q2
        else   
          q(k)=1.e-6   ! a small value for q
        endif
      enddo 
!
      end subroutine pert_rh2q
!
!
      end module m_obs_pert