da_1d_eigendecomposition.inc

1 subroutine da_1d_eigendecomposition( bx, e, l )
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    real, intent(in)         :: bx(:,:)          ! Global vert. background error.
12    
13    real, intent(out)        :: e(:,:)           ! Eigenvectors of Bx.
14    real, intent(out)        :: l(:)             ! Global eigenvalues of Bx.
15    
16    integer                  :: kz               ! Size of 3rd dimension.
17    integer                  :: m                ! Loop counters
18    integer                  :: work             ! Size of work array.
19    integer                  :: info             ! Info code.
20    
21    real, allocatable        :: ecopy(:,:)
22    real, allocatable        :: lcopy(:)
23    real, allocatable        :: work_array(:)
24 
25    if (trace_use_dull) call da_trace_entry("da_1d_eigendecomposition")
26 
27    !-------------------------------------------------------------------------
28    ! [1.0]: Initialise:
29    !-------------------------------------------------------------------------
30    
31    kz = size(bx, dim=1)
32    
33    !-------------------------------------------------------------------------
34    ! [5.0]: Perform global eigenvalue decomposition using LAPACK software:
35    !-------------------------------------------------------------------------
36    
37    work = 3 * kz - 1
38    allocate( work_array(1:work) )
39    
40    allocate( ecopy(1:kz,1:kz) )
41    allocate( lcopy(1:kz) )
42    
43    ecopy(:,:) = bx(:,:)
44 
45    lcopy = 0.0
46 
47    call dsyev( 'V', 'U', kz, ecopy, kz, lcopy, &
48               work_array, work, info )
49    
50    if ( info /= 0 ) then
51       write(unit=message(1),fmt='(A,I4,A)') &
52          ' da_1d_eigendecomposition: info = ', &
53          info,' - error in decomposition.'
54       call da_error(__FILE__,__LINE__,message(1:1))
55    end if
56    
57    !--Swap order of eigenvalues, vectors so 1st is one with most
58    !  variance, etc:
59    
60    do m=1,kz
61       l(m) = lcopy(kz+1-m)
62       e(1:kz,m) = ecopy(1:kz,kz+1-m)
63    end do
64    
65    deallocate (work_array)
66    deallocate (ecopy)
67    deallocate (lcopy)
68 
69    if (trace_use_dull) call da_trace_exit("da_1d_eigendecomposition")
70    
71 end subroutine da_1d_eigendecomposition
72 
73