var/da/da_tools/da_eof_decomposition.inc

   1 subroutine da_eof_decomposition (kz, bx, e, l)
   2
   3    !---------------------------------------------------------------------------
   4    ! Purpose: Compute eigenvectors E and eigenvalues L of vertical covariance
   5    !          matrix
   6    !          B_{x} defined by equation:  E^{T} B_{x} E = L, given input kz x kz
   7    !          BE field.
   8    !---------------------------------------------------------------------------
   9
  10    implicit none
  11
  12    integer, intent(in)  :: kz               ! Dimension of error matrix.
  13    real,    intent(in)  :: bx(1:kz,1:kz)    ! Vert. background error.
  14    real*8,  intent(out) :: e(1:kz,1:kz)     ! Eigenvectors of Bx.
  15    real*8,  intent(out) :: l(1:kz)          ! Eigenvalues of Bx.
  16
  17    integer :: work             ! Size of work array.
  18    integer :: m                ! Loop counters
  19    integer :: info             ! Info code.
  20
  21    real*8  :: work_array(1:3*kz-1)
  22    real*8  :: ecopy(1:kz,1:kz)
  23    real*8  :: lcopy(1:kz)
  24
  25    if (trace_use) call da_trace_entry("da_eof_decomposition")
  26
  27    !-------------------------------------------------------------------------
  28    ! [5.0]: Perform global eigenvalue decomposition using LAPACK software:
  29    !-------------------------------------------------------------------------
  30
  31    work = 3 * kz - 1
  32    ecopy(1:kz,1:kz) = bx(1:kz,1:kz)
  33    lcopy(1:kz) = 0.0
  34
  35    call dsyev( 'V', 'U', kz, ecopy, kz, lcopy, work_array, work, info )
  36
  37    if ( info /= 0 ) then
  38       write(unit=message(1),fmt='(A,I4)') &
  39          "Error in decomposition, info = ", info
  40       call da_error(__FILE__,__LINE__,message(1:1))
  41    end if
  42
  43    ! Swap order of eigenvalues, vectors so 1st is one with most variance:
  44
  45    do m = 1, kz
  46       l(m) = lcopy(kz+1-m)
  47       e(1:kz,m) = ecopy(1:kz,kz+1-m)
  48    end do
  49
  50    if (trace_use) call da_trace_exit("da_eof_decomposition")
  51
  52 end subroutine da_eof_decomposition
  53
  54