      subroutine perconuw2 (n,m,l,lh,hx,dcov,y,hc,lc,dc,lm,pxc,sxc,
     &                      dcom,x,dcon,hxt,p,b1,b2,b3)
c... purpose:
c            consider covariance analysis. to compute perturbation
c              matrix, sensitivity matrix, computed covariance
c              and considered covariance matrices, given coefficient
c              matrices hx,hc, one dimensional data  cov matrix dcov,
c              and covariance matrix of the considered parameters dc.
c...author:
c        xiaoping wu
c
c...date: created in Jan, 2003 from perconuw
c         to add leading dimensions of matrices to subroutine arguments
c...
c...parameters:
c      input:
c        n,m,l  --- number of observation, computed parameters and con-
c                   siderded parameters
c        lh         leading dimension of hx,hc matrices in the calling routine
c        hx(n,m)    measurement  matrix for x
c        dcov(n)    one dimensional data noise cov matrix
c        y(n)       observations
c        dcom(m,m)  a priori estimated parameter covariance matrix
c        hc(n,l)    measurement matrix 
c        lc         leading dimension of matrix dc
c        dc(l,l)    covariance matrix of c
c        lm         leading dimension of dcom,dcon,sxc,pxc matrices
c      output:
c        pxc(m,l)   perturbation matrix
c        sxc(m,l)   sensitivity matrix
c        dcom(m,m)  computed covariance matrix
c        x(m)       parameter estimates
c        dcon(m,m)  considered covariance matrix
c        bi and hxt     --- working units
	implicit double precision (a-h,o-z)
	include 'stdunits.h'
	integer dmax,pmax,pele
	parameter (dmax=8500)
c	parameter (dmax=530000 )
	parameter (pmax=600,pele=(pmax+1)*(pmax+2)/2 )
	dimension hx(lh,m),hc(lh,l),dcov(lh),y(lh),p(lh)
	dimension dcom(lm,m),pcov(pele),dc(lc,l)
	dimension hxt(lm,n)
	dimension b1(lh+1,m+1), b2(lm,l), b3(lc,l)
	dimension pxc(lm,l),sxc(lm,l),x(m),dcon(lm,m)

      character*120 rcsID
      data rcsID
     &/'$Id: perconuw2.f,v 1.8 2005/10/08 02:25:16 xpw Exp $'/

	call TZERO(pcov,m+1,1,m+1)
c. convert a priori parameter covariance into one dimensional array
	icount = 0
	do j=1,m
        do i=1,j
	icount = icount + 1
	pcov(icount) = dcom(i,j)
	end do
	end do

	sum=0.d0
	do j=1,m
	sum=sum+dcom(j,j)
	end do
c. augment the measurement matrix
	do i=1,n
	do j=1,m
	b1(i,j)=hx(i,j)/sqrt(dcov(i))
	end do
	end do
c 
	do j=1,m
	b1(n+1,j)=0.d0
	end do

	do i=1,n
	b1(i,m+1)=y(i)/sqrt(dcov(i))
	end do

	b1(n+1,m+1)=1.d0

	call COV2UD(pcov,m,ierr,stderr)
	call RUDR(pcov,pcov,-m,1,cnb)
c...cov2ud and rudr handles null pcov matrix gracefully
c   here is another protection against null dcom matrix
        if(sum.eq.0.d0) call TZERO(pcov,m+1,1,m+1)
	if(abs(sum).lt.1.d-15 .and. sum.ne.0.d0)then
	write(stderr,*)' non-zero but very small a priori covariance'
	stop
	end if

	rsos = 0.d0
	call THHC(pcov,m,b1,lh+1,n,rsos,1,m)
      WRITE(stdout,'(1X,''CHI^2      =  '',E10.4)') rsos*rsos
      WRITE(stdout,'(1X,''DOF        =  '',I10)') n - m
      nof = n - m
      IF (nof .NE. 0.D0) THEN
       WRITE(stdout,'(1X,''CHI^2/DOF  =  '',F10.2)') rsos*rsos/nof
      END IF
c. convert to covariance
	call RUDR(pcov,pcov,m,1,cnb)
	call UD2COV(pcov,pcov,m)

	icount = 0
        do j= 1,m
         do i =1, j
          icount = icount + 1
          dcom(i,j) = pcov(icount)
          dcom(j,i) = dcom(i,j)
         end do
        end do

        do i = 1, m
         icount = icount + 1
         x(i) = pcov(icount)
        end do
	if(l.eq.0)return

c. get AT P
 	do i=1,n
 	p(i)=1.d0/dcov(i)
 	end do
c
	do i=1,m
	do j=1,n
	hxt(i,j)=hx(j,i)*p(j)
	end do
	end do

c. compute sensitivity matrix
	call ab2cn(m,n,l,lm,lh,lm,hxt,hc,b2)
	call ab2cn(m,m,l,lm,lm,lm,dcom,b2,sxc)
	do 70 i=1,m
	do 70 j=1,l
 70	sxc(i,j)=-sxc(i,j)
c. compute perturbation matrix
	do 50 i=1,l
	do 50 j=1,l
 50	b3(i,j)=0.d0
	do 60 i=1,l
 60	b3(i,i)=dsqrt(dc(i,i))
	call ab2cn(m,l,l,lm,lc,lm,sxc,b3,pxc)
c. compute consider covariance matrix
	call ab2cn(m,l,l,lm,lc,lm,sxc,dc,b2)
	call abt2cn(m,l,m,lm,lm,lm,b2,sxc,dcon)
	do 90 i=1,m
	do 90 j=1,m
 90	dcon(i,j)=dcon(i,j)+dcom(i,j)
	return
	end
c
