!----------------------- EOF分解子程序 -----------------------------------------------!
!要求输入的有:                                                                        !
!   1. m:空间点、网格点或测站;                                                        !
!   2. n:时间(序列长度);                                                              !
!   3. itp:要求输出特征向量和时间系数的个数(itp<=m);                                  !
!   4. ki:对原始场的处理(= -1表示用原始场;= 0表示用距平场;= 1表示用标准化场);         !
!   5. raX(m,n):原始资料矩阵;                                                         !
!输出的有:                                                                            !
!   1. raEgvt(m,m):每一列为一个空间型;存放在file_eof和file_egvt中                     !
!   2. raTime(m,n):每一行为对应空间型的时间系数;存放在file_eof和file_time中;          !
!                  应注意存放在文件时为转置存放,即每列为对应空间型的时间系数;         !
!注意:                                                                                !
!   1. 第33,34行的输出格式与itp的取值有关;                                            !
!   2. 本程序调用了fortran 5.0的smaths.lib中的内部数学函数;                           !
!   3. raA(m,m):协方差矩阵-----动态数组;                                              !
!   4. raEgvl(m,4):第1列为特征值,第2列为累计特征值,                                   !
!                  第3列为方差贡献,第4列为累计方差贡献;-----动态数组(存放在file_eof)  !
!   5. 本程序用到了sub_standard.f90文件;                                              !
!--------------------------------------------------- 程正泉 2000.5 -------------------!
subroutine sub_eof(m,n,itp,ki,raX,               &
                   raEgvt,raTime,                &
                   file_eof,file_egvt,file_time)
implicit none
integer,intent(in)::m,n,itp,ki
real,dimension(m,n),intent(in)::raX
real,dimension(m,m),intent(out)::raEgvt,raTime(m,n)
real,allocatable,dimension(:,:)::raA,raEgvl
real::lamda,temp
integer::i,j
character*50 file_eof,file_egvt,file_time

open(90,file=file_eof);open(91,file=file_egvt);open(92,file=file_time)
101 format(6i15)                         ! 格式随itp的取值而变动
102 format(i4,6f15.4)                    ! 格式随itp的取值而变动

!----- 0.资料处理 ------------------------!
call sub_standard2(m,n,raX,ki)

!----- 1.计算协方差矩阵raA=raX*raX' ------!
allocate(raA(m,m))
call mxytf(m,n,raX,m,              &
             m,n,raX,m,            &
		    m,m,raA,m)
do i=1,m; do j=1,m
   raA(i,j)=raA(i,j)/float(n)
enddo; enddo

!----- 2.计算raA的特征值和特征向量 -------!
allocate(raEgvl(m,4))
call evcsf(m,raA,m,raEgvl(1:m,1),raEgvt,m)
deallocate(raA)

!----- 3.计算时间系数raTime=raEgvt'*raX --!
call mxtyf(m,m,raEgvt,m,           &          
             m,n,raX,m,            &
			   m,n,raTime,m)

!----- 4.验证分解是否正确 ----------------!
write(90,*)'当下面两列值相等或非常近似时,EOF分解正确;否则,分解不正确'
write(90,'(/a8,6x,a12,7x,a22)')'序号','特征值','时间系数平方累加的平均'
do i=1,m
   lamda=0.0
   do j=1,n
      lamda=lamda+raTime(i,j)**2/float(n)
   enddo
   write(90,'(i7,2f20.5)')i,raEgvl(i,1),lamda
enddo
write(90,'(/a43)')'!-----!-----!-----!-----!-----!-----!-----!'

!----- 5.输出EOF分解的结果 ---------------!
write(90,'(a4,8x,a6,6x,a10,4x,a8,4x,a12)')'序号','特征值','累计特征值','方差贡献','累计方差贡献'
temp=0.0
do i=1,m
   temp=temp+raEgvl(i,1)
   raEgvl(i,2)=temp
enddo
do i=1,m
   raEgvl(i,3)=raEgvl(i,1)/raEgvl(m,2)   
enddo
do i=1,m
   raEgvl(i,4)=0.0
   do j=1,i
      raEgvl(i,4)=raEgvl(i,4)+raEgvl(j,3)
   enddo
enddo
do i=1,m
   write(90,'(i3,f15.3,f15.3,f12.3,f14.3)')i,(raEgvl(i,j),j=1,4)
enddo
deallocate(raEgvl)
write(90,'(/a43)')'!-----!-----!-----!-----!-----!-----!-----!'
write(90,*)'空间型:每一列为一个空间型'
write(90,101)(i,i=1,itp)
do i=1,m
   write(90,102)i,(raEgvt(i,j),j=1,itp)
   write(91,102)i,(raEgvt(i,j),j=1,itp)
enddo
write(90,'(/a43)')'!-----!-----!-----!-----!-----!-----!-----!'
write(90,*)'时间系数:每一列为对应空间型的时间系数'
write(90,101)(i,i=1,itp)
do j=1,n
   write(90,102)j,(raTime(i,j),i=1,itp)
   write(92,102)i,(raTime(i,j),i=1,itp)
enddo

end