var/convertor/average_be/da_1d_eigendecomposition.F

   1 subroutine da_1d_eigendecomposition( bx, e, l, k )
   2
   3    !------------------------------------------------------------------------------
   4    !  Purpose: Compute eigenvectors E and eigenvalues L of vertical covariance matrix
   5    !           B_{x} defined by equation:  E^{T} B_{x} E = L, given input 3D field of
   6    !           errors (sum over all horizontal locations).
   7    !------------------------------------------------------------------------------
   8
   9    implicit none
  10
  11    integer, intent(in)               :: k
  12
  13    real, dimension(k,k), intent(in)  :: bx      ! Global vert. background error.
  14    real, dimension(k,k), intent(out) :: e       ! Eigenvectors of Bx.
  15    real, dimension(k), intent(out)   :: l       ! Global eigenvalues of Bx.
  16
  17    integer                  :: kz               ! Size of 3rd dimension.
  18    integer                  :: m                ! Loop counters
  19    integer                  :: work             ! Size of work array.
  20    integer                  :: info             ! Info code.
  21
  22    real, allocatable        :: ecopy(:,:)
  23    real, allocatable        :: lcopy(:)
  24    real, allocatable        :: work_array(:)
  25
  26 !   if (trace_use_dull) call da_trace_entry("da_1d_eigendecomposition")
  27
  28    !-------------------------------------------------------------------------
  29    ! [1.0]: Initialise:
  30    !-------------------------------------------------------------------------
  31
  32    kz = size(bx, dim=1)
  33
  34    !-------------------------------------------------------------------------
  35    ! [5.0]: Perform global eigenvalue decomposition using LAPACK software:
  36    !-------------------------------------------------------------------------
  37
  38    work = 3 * kz - 1
  39    allocate( work_array(1:work) )
  40
  41    allocate( ecopy(1:kz,1:kz) )
  42    allocate( lcopy(1:kz) )
  43
  44    ecopy(:,:) = bx(:,:)
  45
  46    lcopy = 0.0
  47
  48    call dsyev( 'V', 'U', kz, ecopy, kz, lcopy, &
  49               work_array, work, info )
  50
  51 !   if ( info /= 0 ) then
  52 !      write(unit=message(1),fmt='(A,I4,A)') &
  53 !         ' da_1d_eigendecomposition: info = ', &
  54 !         info,' - error in decomposition.'
  55 !      call da_error(__FILE__,__LINE__,message(1:1))
  56 !   end if
  57
  58    !--Swap order of eigenvalues, vectors so 1st is one with most
  59    !  variance, etc:
  60
  61    do m=1,kz
  62       l(m) = lcopy(kz+1-m)
  63       e(1:kz,m) = ecopy(1:kz,kz+1-m)
  64    end do
  65
  66    deallocate (work_array)
  67    deallocate (ecopy)
  68    deallocate (lcopy)
  69
  70 !   if (trace_use_dull) call da_trace_exit("da_1d_eigendecomposition")
  71
  72 end subroutine da_1d_eigendecomposition