subroutine da_varbc_precond (iv)

   !---------------------------------------------------------------------------
   !  PURPOSE:  Calculate covariance matrix b/w bias predictors 
   !            for VarBC preconditioning
   !
   !  Called from da_get_innov_vector_radiance
   !
   !  HISTORY: 11/07/2007 - Creation                     Tom Auligne
   !---------------------------------------------------------------------------

   implicit none

   type (iv_type), intent (inout)   :: iv             ! Innovation

   integer                          :: inst, n, i, j, k, ii, jj
   integer                          :: npred, npredmax, num_rad, num_rad_active
   real, allocatable 		    :: hessian(:,:)
   real*8, allocatable 		    :: eignvec(:,:), eignval(:)
   real                             :: hessian_local, bgerr_local, pred_i, pred_j
   
   if (trace_use) call da_trace_entry("da_varbc_precond")
   
   if ( iv%num_inst < 1 ) RETURN
   
   write(unit=stdout,fmt='(A)') 'VARBC: Estimate Hessian for preconditioning'

   do inst = 1, iv%num_inst      ! loop for sensors
      npredmax = iv%instid(inst)%varbc_info%npredmax
 
      allocate ( hessian(npredmax, npredmax) )
      allocate ( eignvec(npredmax, npredmax) )
      allocate ( eignval(npredmax)           )
	 
      do k = 1, iv%instid(inst)%nchan         ! loop for channels
         npred = iv%instid(inst)%varbc(k)%npred
         if (npred   <= 0) cycle	      ! VarBC channel only     

         num_rad        = iv%instid(inst)%num_rad
         num_rad_active = COUNT( (iv%instid(inst)%info%proc_domain(1,1:num_rad)) &
	                    .AND.(iv%instid(inst)%tb_qc(k,1:num_rad) > qc_varbc_bad) )
         iv%instid(inst)%varbc(k)%nobs  = wrf_dm_sum_integer(num_rad_active)
	    
         if ( satinfo(inst)%iuse(k) == 1 ) &
            write(unit=stdout,fmt='(A,I6,3A,I5)') &
	    'VARBC:',iv%instid(inst)%varbc(k)%nobs,' active observations for ', &
            trim(adjustl(iv%instid(inst)%rttovid_string)),' channel',           &
	    iv%instid(inst)%ichan(k)
	    
         if (iv%instid(inst)%varbc(k)%nobs == 0) cycle
	 
        !---------------------------------------------------------	 
	! Calculate estimation of the Hessian for preconditioning
        !---------------------------------------------------------	 
         do i = 1, npred
	    ii = iv%instid(inst)%varbc(k)%ipred(i)

	   ! Observation term
           !------------------		 
            do j = i, npred
	       jj = iv%instid(inst)%varbc(k)%ipred(j)
	       hessian_local = 0.0	
	       
               do n= 1, num_rad      ! loop for pixel      
                  if (iv%instid(inst)%tb_qc(k,n) <= qc_varbc_bad)   cycle  ! good obs only
	          if (.not. iv%instid(inst)%info%proc_domain(1,n))  cycle  ! do not sum up HALO data
		  
		  if (ii == iv%instid(inst)%varbc_info%gammapred) then
		     pred_i = iv%instid(inst)%gamma_jacobian(k,n)
		  else
   	    	     pred_i = iv%instid(inst)%varbc_info%pred(ii,n)
		  end if   

		  if (jj == iv%instid(inst)%varbc_info%gammapred) then
		     pred_j = iv%instid(inst)%gamma_jacobian(k,n)
		  else
		     pred_j = iv%instid(inst)%varbc_info%pred(jj,n)		  
		  end if   
			      
       	          hessian_local = hessian_local + pred_i * pred_j / &
	                          iv%instid(inst)%tb_error(k,n)**2				     
               end do                                               !  end loop for pixel
	       
	       ! Sum hessian preconditioning across processors
               hessian(i,j) = wrf_dm_sum_real(hessian_local)	
	       hessian(j,i) = hessian(i,j)  	       
	    end do
         
	   ! Background term
           !-----------------
            if (iv%instid(inst)%varbc_info%nbgerr(ii) <= 0) cycle
	    bgerr_local = 0.0
	    do n= 1, num_rad    
               if (iv%instid(inst)%tb_qc(k,n) <= qc_varbc_bad)   cycle  ! good obs only
               if (.not. iv%instid(inst)%info%proc_domain(1,n))  cycle  ! do not sum up HALO data

   	       bgerr_local = bgerr_local + iv%instid(inst)%tb_error(k,n)**2  / &
                                             varbc_nbgerr
	                                   ! iv%instid(inst)%varbc_info%nbgerr(ii)
	    end do

            iv%instid(inst)%varbc(k)%bgerr(i) = wrf_dm_sum_real(bgerr_local) / &
	                                        iv%instid(inst)%varbc(k)%nobs

	    if (iv%instid(inst)%varbc(k)%bgerr(i) > 0) &          
   	       hessian(i,i) = hessian(i,i) + 1/iv%instid(inst)%varbc(k)%bgerr(i)
	 end do   

        !--------------------------------------------------	 
	! Preconditioning = inverse square root of Hessian
        !--------------------------------------------------	 	 
         hessian = hessian / varbc_factor**2
	 
	 if (npred == 1) then
	    iv%instid(inst)%varbc(k)%vtox(1,1)     = 1.0 / sqrt(hessian(1,1))
         else 
	    call da_eof_decomposition(npred, hessian(1:npred,1:npred), &
	    	            eignvec(1:npred,1:npred),eignval(1:npred))

	    if (ANY(eignval(1:npred) <= 0)) then		    
	       write(unit=stdout,fmt='(3A,I4,A,10F18.5)') &
	       'VARBC: non-positive Hessian for ', iv%instid(inst)%rttovid_string, &
	              ' ,channel ',k,'--> Eigenvalues =',eignval(1:npred) 
	       do i = 1, npred
	          if (hessian(i,i) > 0) &
	             iv%instid(inst)%varbc(k)%vtox(i,i) = 1.0 / sqrt(hessian(i,i))
	       end do
	    else
 	       do i = 1, npred
	          do j = i, npred
                     iv%instid(inst)%varbc(k)%vtox(i,j) = sum(          &
		                           eignvec(i,1:npred)         * &
			                   sqrt(1.0/eignval(1:npred)) * &
			                   eignvec(j,1:npred) )
		       
		     iv%instid(inst)%varbc(k)%vtox(j,i) = &
          	     iv%instid(inst)%varbc(k)%vtox(i,j)
                  end do
	       end do
	    end if
	 end if   
!         do i = 1, npred
!	    print*,inst,k,i,(iv%instid(inst)%varbc(k)%vtox(i,j),j=1,npred)
!	 end do
      end do                     !  end loop for channels	
      deallocate(hessian, eignvec, eignval)   
   end do                        !  end loop for sensor


   if (trace_use) call da_trace_exit("da_varbc_precond")

 end subroutine da_varbc_precond