s_lin_eof Subroutine

subroutine s_lin_eof(x, nd, nv, pc, eof, ew, opt, wt, r2, eof_scaled)

! ==== Description
!! Empirical Orthogonal Function (EOF) analysis

! ==== Declarations
  integer(i4), intent(in)            :: nd                !! number of rows
  integer(i4), intent(in)            :: nv                !! number of columns
  real(wp)   , intent(in)            :: x(nd,nv)          !! input data
  integer(i4), intent(in) , optional :: opt               !! 0 = covariance, 1 = correlation
  real(wp)   , intent(in) , optional :: wt(nv)            !! optional weights (default = 1.0/n)
  real(wp)   , intent(out)           :: pc(nd,nv)         !! principal components
  real(wp)   , intent(out)           :: eof(nv,nv)        !! EOFs/eigenvectors (unweighted)
  real(wp)   , intent(out)           :: ew(nv)            !! eigenvalues
  real(wp)   , intent(out), optional :: r2(nv)            !! explained variance (fraction)
  real(wp)   , intent(out), optional :: eof_scaled(nv,nv) !! EOFs/eigenvectors scaled for plotting
  real(wp)                           :: eof_tmp(nv,nv)    !! temp storage for EOFs
  real(wp)                           :: x_w(nd,nv)        !! working copy of data
  integer(i4)                        :: opt_w             !! final value for opt
  real(wp)                           :: wt_w(nv)          !! final values for wt
  real(wp)                           :: c(nv,nv)          !! covariance/correlation matrix
  real(wp)                           :: w(nv)             !! matrix for eigenvalues (diagonal)
  integer(i4)                        :: nn                !! number of nonzero eigenvalues
  integer(i4)                        :: i, j              !! loop indices
  real(wp)                           :: tmp               !! temporary real

! ==== Instructions

  ! ---- handle input

  ! check matrix dimensions
  if (nd .lt. 2 .or. nv .lt. 1) then
     ! write error message and stop if invalid
     call s_err_print(fsml_error(3))
     error stop
  endif

  ! weight options; will default to 1/nv if not specified
  tmp = real(nv, kind=wp)
  wt_w = 1.0_wp / tmp
  if (present(wt)) wt_w = wt

  ! matrix options: covariance (0, default) or correlation
  opt_w = 0
  if (present(opt)) opt_w = opt

  ! make working copy of data
  x_w = x

  ! ---- construct covariance or correlation matrix

  ! prepare data
  do j = 1, nv
     ! centre (get anomalies)
     tmp = f_sts_mean_core(x_w(:,j))
     x_w(:,j) = x_w(:,j) - tmp
     ! standardise if specified (if correlation matrix is to be used)
     if (opt_w .eq. 1) then
        tmp = f_sts_std_core(x_w(:,j), 0.0_wp)
        if (tmp .gt. 0.0_wp) then
           x_w(:,j) = x_w(:,j) * (1.0_wp / tmp)
        else
           call s_err_print("[fsml error] Standard deviation&
                           & is zero for a column.")
           error stop
        endif
     endif
     ! apply weights
     x_w(:,j) = x_w(:,j) * sqrt(wt_w(j))
  enddo

  ! construct matrix
  tmp = real(nd - 1, kind=wp)
  c = matmul(transpose(x_w), x_w) / tmp

  ! ---- calculate outputs

  ! eigen-decomposition using stdlib eigh
  call eigh(c, w, vectors=eof_tmp)

  ! extract and reorder eigenvalues/eigenvectors in descending order
  eof = 0.0_wp
  ew  = 0.0_wp
  nn  = 0
  do i = nv, 1, -1
     if (w(i) .gt. 0.0_wp) then
        nn = nn + 1
        ew(nn) = w(i)
        eof(:,nn) = eof_tmp(:,i)
     endif
  enddo

  ! compute principal components
  pc = 0.0_wp
  pc(:,1:nn) = matmul(x_w, eof(:,1:nn))

  ! undo weight scaling for EOFs
  if (nn .gt. 0) then
     do i = 1, nv
        tmp = 1.0_wp / sqrt(wt_w(i))
        eof(i,1:nn) = eof(i,1:nn) * tmp
     enddo
  endif

  ! scale EOFs for plotting
  if (present(eof_scaled)) then
     if (nn .gt. 0) then
        eof_scaled = 0.0_wp
        do j = 1, nn
           tmp = sqrt(ew(j))
           eof_scaled(:,j) = eof(:,j) * tmp
        enddo
     endif
  endif

  ! explained variance (fraction)
  if (present(r2)) then
     r2 = 0.0_wp
     r2(1:nn) = ew(1:nn) / sum(ew(1:nn))
  endif

end subroutine s_lin_eof