tcorr.f90

求两个序列的相关关系的fortan程序 看看

!尺度分离加 点面相关


program main
PARAMETER (Nt=54,nx=192,ny=94,ki=11)
REAL skt(nx,ny,nt),corr(nx,ny),sm(55),smi(nt)
real v(nt),y(nt),x(nt),t(nt),x1(nt)
real nskt(nx,ny,nt),nsmi(nt),ndskt(nx,ny,nt),ndsmi(nt)
real m(nx,ny,nt),d(nt),mqu(nx,ny,nt),dqu(nt)
character*50 infile1,infile2,outfile


infile1='i:\monsoonindex\wmi.grd'
infile2='i:\skt\skt.grd'
outfile='i:\monsoonindex\corr.grd'

OPEN(11,FILE=infile1,FORM='BINARY')
     READ(11)(sm(i),I=1,55)
CLOSE(11)
write(*,*)(sm(i),i=1,55)
do i=1,nt
smi(i)=sm(i+1)
enddo
write(*,*)(smi(i),I=1,nt)
OPEN(12,FILE=infile2,FORM='BINARY')
     READ(12)(((skt(ix,iy,it),ix=1,nx),iy=1,ny),it=1,nt)
CLOSE(12)
!-------------------chidu fenli
do  ix=1,nx
  do  iy=1,ny
    do it=1,nt
       v(it)=skt(ix,iy,it)
    enddo
    call sub_ph2(nt,ki,v,y) 
    do it=1,nt
	   ndskt(ix,iy,it)=y(it)
       nskt(ix,iy,it)=skt(ix,iy,it)-y(it)
    enddo
  enddo
enddo


!-----------------------------------------       
	  call sub_ph2(nt,ki,smi,ndsmi) 
        do i=1,nt
	    nsmi(i)=smi(i)-ndsmi(i)
	  enddo
!-----------------------------------------
write(*,*) 'nsmi',nsmi
write(*,*) 'ndsmi',ndsmi
write(*,*) 'smi',smi

!===============去趋势后相关====================

!---------------------------------------------
   do ix=1,nx
      do  iy=1,ny
	    do it=1,nt
          m(ix,iy,it)=nskt(ix,iy,it)
        enddo
      enddo
   enddo
  do it=1,nt
     d(it)=nsmi(it)
  enddo
  
 !---------------jisuan xiangguan-----------------------
    do ix=1,nx
	  do  iy=1,ny
	    do it=1,nt
	    x(it)=m(ix,iy,it)
	    t(it)=real(it+1947)
	    enddo
	    call xianxing(x,t,nt,a,b,r)
	    do it=1,nt
	    x1(it)=a+b*t(it)
	    mqu(ix,iy,it)=x(it)-x1(it)
        enddo
	  enddo
   enddo
      
	    do it=1,nt
	    x(it)=d(it)
	    t(it)=real(it+1948)
	    enddo
	    call xianxing(x,t,nt,a,b,r)
	    do it=1,nt
	    x1(it)=a+b*t(it)
	    dqu(it)=x(it)-x1(it)
        enddo

	do ix=1,nx 
	  do iy=1,ny
	    do it=1,nt
	      x(it)=mqu(ix,iy,it)
	    enddo
        call corrr(nt,dqu,x,r)
	    corr(ix,iy)=r
	  enddo
	enddo
!-----------------------------------------
	write(*,*)'corr=', ((corr(ix,iy),ix=1,nx),iy=1,ny)
!-----------------------------------
open(10,file=outfile,form='binary')
write(10)((corr(ix,iy),ix=1,nx),iy=1,ny)
close(10)
end



!!!==============zichengxu================

	  
!=================子程序1：相关 ====================================      
	subroutine corrr(n,x,y,r)                                     
      dimension x(n),y(n)
	ax=0.0
	ay=0.0
	do i=1,n
	ax=ax+x(i)
	ay=ay+y(i)   
	enddo
	ax=ax/real(n)
	ay=ay/real(n)

      sxy=0.0
	sx2=0.0
	sy2=0.0
	do i=1,n
	sxy=sxy+(x(i)-ax)*(y(i)-ay)  
	sx2=sx2+(x(i)-ax)*(x(i)-ax)
	sy2=sy2+(y(i)-ay)*(y(i)-ay) 
	enddo
	r=sxy/(sqrt(sx2)*sqrt(sy2))
	end        
      
!----------------- 平滑子程序 --------------------------- !
!本文件共有两个子程序:                                    !
!       1.sub_ph1(n,ki,x,y):线性平滑;                     !
!       2.sub_ph2(n,ki,x,y):二项式平滑;                   !
!要求输入的有:                                            !
!       1.x(n):长度为n的序列;                             !
!       2.ki:平滑步长;                                    !
!输出的有:                                                !
!       1.y(n):平滑后的序列;                              !
!------------------------------------- 程正泉 2000.5 -----!
!=========================================================!

!=========================================================!
subroutine sub_ph2(n,ki,x,y)                              !
implicit none                                             !
integer,intent(in)::n,ki                                  !
real,dimension(n),intent(in)::x                           !
real,dimension(n),intent(out)::y                          !
real,allocatable,dimension(:)::h                          !
integer::i,j                                              !
real::temp                                                !
                                                          !
allocate(h(-ki/2:ki/2))                                   !
call coef(h,ki)                                           !
do i=1,n                                                  !
  if(i<=ki/2)then                                         !
    y(i)=0.0; temp=0.0                                    !
    do j=1-i,ki/2                                         !
       y(i)=y(i)+h(j)*x(i+j)                              !
       temp=temp+h(j)                                     !
    enddo                                                 !
    y(i)=y(i)/temp                                        !
  elseif((i > ki/2).and.(i < n-Ki/2+1))then               !
    y(i)=0.0; temp=0.0                                    !
    do j=-ki/2,ki/2                                       !
       y(i)=y(i)+h(j)*x(i+j)                              !
       temp=temp+h(j)                                     !
    enddo                                                 !
    y(i)=y(i)/temp                                        !
  else                                                    !
    y(i)=0.0; temp=0.0                                    !
    do j=-ki/2,n-i                                        !
       y(i)=y(i)+h(j)*x(i+j)                              !
       temp=temp+h(j)                                     !
    enddo                                                 !
    y(i)=y(i)/temp                                        !
  endif                                                   !
enddo                                                     !
deallocate(h)                                             !
end                                                       !
!     ---     ---     ---     ---    ---    ---    ---    !
 subroutine coef(h,k)                                      !
!二项式的系数,参考丁裕国教授所编著的气象...               !
      implicit none                                             !
      integer,intent(in)::k                                     !
      real,dimension(-k/2:k/2),intent(out)::h                   !
      integer,external::ppp                                     !
      integer::i                                                !
      do i=-k/2,k/2                                             !
      h(i)=ppp(k-1)*1.0/(ppp(i+k/2)*ppp(k-1-i-k/2))           !
      enddo                                                     !
      end                                                       !
!     ---     ---     ---    ---     ---     ---    ---   !
!求n的阶乘(n!)                                            !
      integer function ppp(n)                                   !
      implicit none                                             !
      integer,intent(in)::n                                     !
      integer::i                                                !
      if(n == 0)then                                            !
      ppp=1                                                   !
      else                                                      !
      ppp=1                                                   !
      do i=1,n                                                !
      ppp=ppp*i                                             !
      enddo                                                   !
      endif                                                     !
      end                                                       !
!=========================================================!


!===================zichengxu3:xianxing qushi==================================	
	subroutine xianxing(x,t,nt,a,b,r)
      real x(nt),t(nt),a,b,r

	sxt=0.0
	sx=0.0
	st=0.0
	stt=0.0
	do it=1,nt
      sxt=x(it)*t(it)+sxt
	st=t(it)+st
	sx=x(it)+sx
	stt=t(it)*t(it)+stt
	enddo
	xa=sx/real(nt)
	ta=st/real(nt)

	sax=0.0
	sat=0.0
	saxt=0.0
	do it=1,nt
	sax=sax+(x(it)-xa)*(x(it)-xa)
	sat=sat+(t(it)-ta)*(t(it)-ta)
	saxt=saxt+(x(it)-xa)*(t(it)-ta)
	enddo
       r=saxt/(sqrt(sax*sat))  	    
	 b=(sxt-sx*st/real(nt))/(stt-st*st/real(nt))
	 a=sx/real(nt)-b*st/real(nt)
      end