BEOF.f90

!c-----*----------------------------------------------------6---------7--
!C     EMPIRICAL ORTHOGONAL FUNCTIONS (EOF's)
!c     This subroutine applies the EOF approach to analysis time series 
!c       of meteorological field f(m,n).
!c     input: m,n,mnl,f(m,n),ks
!c       m: number of grid-points
!c       n: lenth of time series
!c       mnl=min(m,n)
!c       f(m,n): the raw spatial-temporal seires
!c       ks: contral parameter
!c           ks=-1: self, i.e., for the raw time series; 
!c           ks=0: departure, i.e., for the anomaly time series from climatological mean;
!c           ks=1: normalized departure, i.e., for the normalized time series; 
!c     output: egvt,ecof,er
!c       egvt(m,mnl): array of eigenvactors
!c       ecof(mnl,n): array of time coefficients for the respective eigenvectors
!c       er(mnl,1): lamda (eigenvalues), its sequence is from big to small.
!c       er(mnl,2): accumulated eigenvalues from big to small
!c       er(mnl,3): explained variances (lamda/total explain) from big to small
!c       er(mnl,4): accumulated explaned variances from big to small
!c     work variables:

!c     Last updated by Dr. Jianping Li on October 20, 2005.

!c     Citation: Li, J., 2005: EOF subroutine, http://web.lasg.ac.cn/staff/ljp/subroutine/EOF.f.

!c-----
    subroutine eof(m,cov,Eval)
        integer::m
        real:: er(m,4)
        real:: cov(m,m),Eval(m),s(m,m),d(m) !work array
!c---- Eigenvalues and eigenvectors by the Jacobi method 
      call jacobi(m,cov,s,d,0.00001)
!c---- Specified eigenvalues 
      call arrang(m,d,s,er)
      cov=s
      Eval(:)=er(:,1)
    return
    end subroutine
!c-----*----------------------------------------------------6---------7--
!c     Preprocessing data to provide a field by ks.
!c     input: m,n,f
!c       m: number of grid-points
!c       n: lenth of time series
!c       f(m,n): the raw spatial-temporal seires
!c       ks: contral parameter
!c           ks=-1: self, i.e., for the raw time series; 
!c           ks=0: departure, i.e., for the anomaly time series from climatological mean;
!c           ks=1: normalized departure, i.e., for the normalized time series; 
!c     output: f
!c       f(m,n): output field based on the control parameter ks.
!c     work variables: fw(n)
    subroutine transf(m,n,f,ks)
      real:: f(m,n)
      real:: fw(n),wn(m)           !work array
      real::af,sf,vf
	i0=0
	do i=1,m
	  do j=1,n
          fw(j)=f(i,j)
        enddo
	  call meanvar(n,fw,af,sf,vf)
	  if(sf.eq.0.)then
	    i0=i0+1
	    wn(i0)=i
	  endif
	enddo
	if(i0.ne.0)then
	  write(*,*)'****  FAULT  ****'
	  write(*,*)' The program cannot go on because& 
          the original field has invalid data.'
	  write(*,*)' There are totally ',i0,&
          'gridpionts with invalid data.'
     	  write(*,*)' These invalid data are those whose standard variance& 
          equal zero.'
     	  write(*,*)' The array WN stores the positions of those invalid& 
          grid-points. You must pick off those invalid data from the orignal&
          field and then reinput a new field to calculate its EOFs.'   
	  write(*,*)'****  FAULT  ****'
	  stop
	endif	    
      if(ks.eq.-1)return
      if(ks.eq.0)then                !anomaly of f
        do i=1,m
          do j=1,n
            fw(j)=f(i,j)
          enddo
          call meanvar(n,fw,af,sf,vf)
          do j=1,n
            f(i,j)=f(i,j)-af
          enddo
        enddo
        return
      endif
      if(ks.eq.1)then                 !normalizing f
        do i=1,m
          do j=1,n
            fw(j)=f(i,j)
          enddo
          call meanvar(n,fw,af,sf,vf)
          do j=1,n
            f(i,j)=(f(i,j)-af)/sf
          enddo
        enddo
      endif
      return
      end subroutine
!c-----*----------------------------------------------------6---------7--
!c     Crossed product martix of the data. It is n times of 
!c       covariance matrix of the data if ks=0 (i.e. for anomaly). 
!c     input: m,n,mnl,f
!c       m: number of grid-points
!c       n: lenth of time series
!c       mnl=min(m,n)
!c       f(m,n): the raw spatial-temporal seires
!c     output: cov(mnl,mnl)  
!c       cov(m,n)=f*f' or f'f dependes on m and n.
!c         It is a mnl*mnl real symmetric matrix.
      subroutine crossproduct(m,n,mnl,f,cov)
      real:: f(m,n),cov(mnl,mnl)
      if(n-m) 10,50,50
  10  do 20 i=1,mnl
      do 20 j=i,mnl
        cov(i,j)=0.0
        do is=1,m
          cov(i,j)=cov(i,j)+f(is,i)*f(is,j)
        enddo
        cov(j,i)=cov(i,j)
  20  continue
      return
  50  do 60 i=1,mnl
      do 60 j=i,mnl
        cov(i,j)=0.0
        do js=1,n
          cov(i,j)=cov(i,j)+f(i,js)*f(j,js)
        enddo
        cov(j,i)=cov(i,j)
  60  continue
      return
      end subroutine
!c-----*----------------------------------------------------6---------7--
!c     Computing eigenvalues and eigenvectors of a real symmetric matrix
!c       a(m,m) by the Jacobi method.
!c     input: m,a,s,d,eps 
!c       m: order of matrix
!c       a(m,m): the covariance matrix
!c       eps: given precision
!c     output: s,d 
!c       s(m,m): eigenvectors
!c       d(m): eigenvalues
      subroutine jacobi(m,a,s,d,eps)
      real:: a(m,m),s(m,m),d(m)
      do 30 i=1,m
      do 30 j=1,i
        if(i-j) 20,10,20
  10    s(i,j)=1.
        go to 30
  20    s(i,j)=0.
        s(j,i)=0.
  30  continue
      g=0.
      do 40 i=2,m
        i1=i-1
        do 40 j=1,i1
  40      g=g+2.*a(i,j)*a(i,j)
      s1=sqrt(g)
      s2=eps/float(m)*s1
      s3=s1
      l=0
  50  s3=s3/float(m)
  60  do 130 iq=2,m
        iq1=iq-1
        do 130 ip=1,iq1
        if(abs(a(ip,iq)).lt.s3) goto 130
        l=1
        v1=a(ip,ip)
        v2=a(ip,iq)
        v3=a(iq,iq)
        u=0.5*(v1-v3)
        if(u.eq.0.0) g=1.
        if(abs(u).ge.1e-10) g=-sign(1.,u)*v2/sqrt(v2*v2+u*u)
        st=g/sqrt(2.*(1.+sqrt(1.-g*g)))
        ct=sqrt(1.-st*st)
        do 110 i=1,m
          g=a(i,ip)*ct-a(i,iq)*st
          a(i,iq)=a(i,ip)*st+a(i,iq)*ct
          a(i,ip)=g
          g=s(i,ip)*ct-s(i,iq)*st
          s(i,iq)=s(i,ip)*st+s(i,iq)*ct
  110     s(i,ip)=g
        do 120 i=1,m
          a(ip,i)=a(i,ip)
  120     a(iq,i)=a(i,iq)
        g=2.*v2*st*ct
        a(ip,ip)=v1*ct*ct+v3*st*st-g
        a(iq,iq)=v1*st*st+v3*ct*ct+g
        a(ip,iq)=(v1-v3)*st*ct+v2*(ct*ct-st*st)
        a(iq,ip)=a(ip,iq)
  130 continue
      if(l-1) 150,140,150
  140 l=0
      go to 60
  150 if(s3.gt.s2) goto 50
      do 160 i=1,m
        d(i)=a(i,i)
  160 continue
      return
      end subroutine
!c-----*----------------------------------------------------6---------7--
!c     Provides a series of eigenvalues from maximuim to minimuim.
!c     input: mnl,d,s
!c       d(mnl): eigenvalues 
!c       s(mnl,mnl): eigenvectors
!c     output: er
!c       er(mnl,1): lamda (eigenvalues), its equence is from big to small.
!c       er(mnl,2): accumulated eigenvalues from big to small
!c       er(mnl,3): explained variances (lamda/total explain) from big to small
!c       er(mnl,4): accumulated explaned variances from big to small
      subroutine arrang(mnl,d,s,er)
      real:: d(mnl),s(mnl,mnl),er(mnl,4)
      tr=0.0
      do 10 i=1,mnl
        tr=tr+d(i)
        er(i,1)=d(i)
  10  continue
      mnl1=mnl-1
      do 20 k1=mnl1,1,-1
      do 20 k2=k1,mnl1
        if(er(k2,1).lt.er(k2+1,1)) then
          c=er(k2+1,1)
          er(k2+1,1)=er(k2,1)
          er(k2,1)=c
          do 15 i=1,mnl
            c=s(i,k2+1)
            s(i,k2+1)=s(i,k2)
            s(i,k2)=c
  15      continue
        endif
  20  continue
      er(1,2)=er(1,1)
      do 30 i=2,mnl
        er(i,2)=er(i-1,2)+er(i,1)
  30  continue
      do 40 i=1,mnl
        er(i,3)=er(i,1)/tr
        er(i,4)=er(i,2)/tr
  40  continue
      return
      end subroutine
!c-----*----------------------------------------------------6---------7--
!c     Provides standard eigenvectors and their time coefficients
!c     input: m,n,mnl,f,s,er
!c       m: number of grid-points
!c       n: lenth of time series
!c       mnl=min(m,n)
!c       f(m,n): the raw spatial-temporal seires
!c       s(mnl,mnl): eigenvectors
!c       er(mnl,1): lamda (eigenvalues), its equence is from big to small.
!c       er(mnl,2): accumulated eigenvalues from big to small
!c       er(mnl,3): explained variances (lamda/total explain) from big to small
!c       er(mnl,4): accumulated explaned variances from big to small
!c     output: egvt,ecof
!c       egvt(m,mnl): normalized eigenvectors
!c       ecof(mnl,n): time coefficients of egvt 
      subroutine tcoeff(m,n,mnl,f,s,er,egvt,ecof)
      real:: f(m,n),s(mnl,mnl),er(mnl,4),egvt(m,mnl),ecof(mnl,n)
      real:: v(mnl)  !work array
      do j=1,mnl
        do i=1,m
          egvt(i,j)=0.
        enddo
        do i=1,n
          ecof(j,i)=0.
        enddo
      enddo
!c-----Normalizing the input eignvectors s
      do 10 j=1,mnl
        c=0.
        do i=1,mnl
          c=c+s(i,j)*s(i,j)
        enddo
        c=sqrt(c)
        do i=1,mnl
          s(i,j)=s(i,j)/c
        enddo
  10  continue
!c-----
      if(m.le.n) then
        do js=1,mnl
        do i=1,m
          egvt(i,js)=s(i,js)
        enddo
        enddo
        do 30 j=1,n
          do i=1,m
            v(i)=f(i,j)
          enddo
          do is=1,mnl
          do i=1,m
            ecof(is,j)=ecof(is,j)+v(i)*s(i,is)
          enddo
          enddo
  30    continue
      else
        do 40 i=1,m
          do j=1,n
            v(j)=f(i,j)
          enddo
          do js=1,mnl
          do j=1,n
            egvt(i,js)=egvt(i,js)+v(j)*s(j,js)
          enddo
          enddo
  40    continue
        do 50 js=1,mnl
          do j=1,n
            ecof(js,j)=s(j,js)*sqrt(abs(er(js,1)))
          enddo
          do i=1,m
            egvt(i,js)=egvt(i,js)/sqrt(abs(er(js,1)))
          enddo
  50    continue
      endif
      return
      end subroutine