var/da/da_tools/da_eof_decomposition_test.inc

   1 subroutine da_eof_decomposition_test (kz, bx, e, l)
   2
   3    !------------------------------------------------------------------------------
   4    ! Purpose:
   5    ! [1] Print eigenvalues:
   6    ! [2] Test orthogonality of eigenvectors - sum_k (e_m(k) e_n(k)) = delta_mn:
   7    ! [3] Test eigenvectors completeness - sum_m (e_m(k1) e_m(k2)) = delta_k1k2:
   8    ! [4] Check B correctness: B = E*L*E^T
   9    !------------------------------------------------------------------------------
  10
  11    implicit none
  12
  13    integer, intent(in) :: kz               ! Dimension of BE matrix
  14    real,    intent(in) :: bx(1:kz,1:kz)    ! Global vert. background error.
  15    real*8,  intent(in) :: e(1:kz,1:kz)     ! Eigenvectors of Bx.
  16    real*8,  intent(in) :: l(1:kz)          ! Eigenvalues of Bx.
  17
  18    integer                  :: k, k1, k2, m     ! Loop counters
  19    real                     :: tot_variance     ! Total variance.
  20    real                     :: cum_variance     ! Cumulative variance.
  21    real                     :: max_off_diag     ! Maximum off-diagonal.
  22
  23    real                     :: work(1:kz,1:kz)  ! 2D Work matrix.
  24    real                     :: bc(1:kz,1:kz)    ! 2D Work matrix.
  25    logical                  :: array_mask(1:kz) ! Array mask for MAXVAL.
  26
  27    if (trace_use) call da_trace_entry("da_eof_decomposition_test")
  28
  29    !-------------------------------------------------------------------------
  30    ! [1] Print eigenvalues:
  31    !-------------------------------------------------------------------------
  32
  33    tot_variance = sum(l(1:kz))
  34    cum_variance = 0.0
  35
  36    write(unit=stdout,fmt='(A)')'  Mode    Eigenvalue     Cumulative Variance      e(k,k)'
  37
  38    do k = 1, kz
  39       cum_variance = cum_variance + l(k)
  40       write(unit=stdout,fmt='(I4,4x,e12.4,10x,f8.4,4x,e12.4)') &
  41             k, l(k), cum_variance / tot_variance, e(k,k)
  42    end do
  43
  44    write(unit=stdout,fmt=*)
  45
  46    call da_array_print( 1, e, 'Global Eigenvectors' )
  47
  48    !-------------------------------------------------------------------------
  49    ! [2] Test orthogonality of eigenvectors - sum_k (e_m(k) e_n(k)) = delta_mn:
  50    !-------------------------------------------------------------------------
  51
  52    write(unit=stdout,fmt='(A)')' Eigenvector orthogonality check:'
  53    write(unit=stdout,fmt='(A)')' Mode     Diagonal         Maximum off-diagonal'
  54
  55    do k1 = 1, kz
  56       do k2 = 1, kz
  57          work(k1,k2) = sum(e(1:kz,k1) * e(1:kz,k2))
  58       end do
  59
  60       array_mask(1:kz) =.true.
  61       array_mask(k1) = .false.
  62       max_off_diag = maxval(abs(work(k1,:)),mask=array_mask(:))
  63       write(unit=stdout,fmt='(I4,4x,1pe12.4,10x,1pe12.4)')k1, work(k1,k1), max_off_diag
  64    end do
  65    write(unit=stdout,fmt=*)
  66
  67    !-------------------------------------------------------------------------
  68    ! [3] Test eigenvectors completeness - sum_m (e_m(k1) e_m(k2)) = delta_k1k2:
  69    !-------------------------------------------------------------------------
  70
  71    write(unit=stdout,fmt='(A)')' Eigenvector completeness check:'
  72    write(unit=stdout,fmt='(A)')' Level    Diagonal         Maximum off-diagonal'
  73
  74    do k1 = 1, kz
  75       do k2 = 1, kz
  76          work(k1,k2) = sum(e(k1,1:kz) * e(k2,1:kz))
  77       end do
  78
  79       array_mask(1:kz) =.true.
  80       array_mask(k1) = .false.
  81       max_off_diag = maxval(abs(work(k1,:)),mask=array_mask(:))
  82       write(unit=stdout,fmt='(I4,4x,1pe12.4,10x,1pe12.4)')k1, work(k1,k1), max_off_diag
  83    end do
  84    write(unit=stdout,fmt=*)
  85
  86    !-------------------------------------------------------------------------
  87    ! [4]  check B correctness: B = E*L*E^T
  88    !-------------------------------------------------------------------------
  89
  90    write(unit=stdout,fmt='(a/a)') &
  91         'real and Calculated B (diagonal)', &
  92         'lvl                 real-B                    Calculated-B'
  93
  94    do k=1,kz
  95       do m=1,kz
  96          work(k,m)=l(k)*e(m,k)
  97          bc(k,m)=0.0
  98       end do
  99    end do
 100
 101    do k1=1,kz
 102       do k2=1,kz
 103          do m=1,kz
 104             bc(k1,k2)=bc(k1,k2)+e(k1,m)*work(m,k2)
 105          end do
 106       end do
 107
 108       write(unit=stdout,fmt='(I5,2F20.5)') k1, bx(k1,k1), bc(k1,k1)
 109    end do
 110
 111    do k2=1,kz
 112       write(unit=stdout, fmt='(a,i4/a)') &
 113            'real and Calculated B (off diagonal):', k2, &
 114            'lvl                 real-B                    Calculated-B'
 115
 116       do k1=1,kz
 117         write(unit=stdout,fmt='(I5,2F20.5)') k1, bx(k1,k2), bc(k1,k2)
 118       end do
 119    end do
 120
 121    if (trace_use) call da_trace_exit("da_eof_decomposition_test")
 122
 123 end subroutine da_eof_decomposition_test
 124
 125