矢量合成分析

非我犹存

想问一下大家，李建平老师的矢量合成分析程序里要怎么读写数据呢，主要是根据指数来确定显著区域。很纠结，子程序如下：
!c-----*----------------------------------------------------6---------7--
!c     For a time ve！!ctor series f(n,2), ！!computing:
!c       (1) mean ve!ctor in the high index years of x(n)
!c       (2) mean ve!ctor in the low index years of x(n)
!c       (3) !composite differen!ce between the mean ve!ctor of f(n,2) in high
!c                     index years of x(n) and its !climatology average
!c       (4) !composite differen!ce between the mean ve!ctor of f(n,2) in low
!c                     index years of x(n) and its !climate average
!c       (5) !composite differen!ce between the mean ve!ctors of f(n,2) in high
!c                     and low index years of x(n)
!c     input: n,x(n),f(n,2),!coefh,!coefl
!c       n: the length of time series
!c       x: !control variable (index)
!c       f: given series
!c       !coefh: !control parameter for high index (i.e., high index years are
!c              those in whi!ch x(i) > !coefh)
!c       !coefl: !control parameter for low index (i.e., low index years are
!c              those in whi!ch x(i) < !coefl)
!c     output: fh(2),fl(2),dh(2),dl(2),dhl(2),tn(5)
!c       fh: the mean ve!ctor of f in high index years of x(n)
!c       fl: the mean ve!ctor of f in low index years of x(n)
!c       dh: !composite differen!ce between the mean ve!ctor of f in high index years of x(n)
!c            and its !climate mean (i.e., high index years minus !cliamte mean).
!c       dl: !composite differen!ce between the mean ve!ctor of f in low index years of x(n)
!c            and its !climate mean (i.e., low index years minus !cliamte mean).
!c       dhl: !composite differen!ce between the mean ve!ctors of f in high and low index years
!c             of x(n) (i.e., high minus low index years)
!c       tn(i,j): tn only equals 1. or 0. !corresponding to signifi!cant differen!ce or not.
!c           tn=1. indi!cates that the differen!ce is signifi!cant.
!c           tn=0. indi!cates that the differen!ce is not signifi!cant.
!c           tn(1,j)~tn(5,j) are !corresponding to the 90%,95%,98%,99% and 99.9% !confident levels.
!c           j=1: siginifi!cant test for dh
!c           j=2: siginifi!cant test for dl
!c           j=3: siginifi!cant test for dhl
!c     Mar!ch 10, 2004 by Jianping Li.
      subroutine differhl1V(n,x,f,coefh,coefl,fh,fl,dh,dl,dhl,tn)
      dimension x(n),f(n,2)
      dimension fh(2),fl(2),dh(2),dl(2),dhl(2),tn(5,3)
      real avef(2),hn(n,2),ln(n,2) !work array
      do j=1,2
        call meanvar(n,f(1,j),avef(j),sf,vf)
      enddo
      do j=1,2
        nh=0
        nl=0
        do k=1,n
          if(x(k).gt.coefh)then
            nh=nh+1
               hn(nh,j)=f(k,j)
          endif
          if(x(k).lt.coefl)then
            nl=nl+1
            ln(nl,j)=f(k,j)
          endif
        enddo
        call meanvar1(nh,hn(1,j),fh(j))
        call meanvar1(nl,ln(1,j),fl(j))
        dh(j)=fh(j)-avef(j)
        dl(j)=fl(j)-avef(j)
        dhl(j)=fh(j)-fl(j)
      enddo
      call diff_t_testV1(n,nh,f(1,1),f(1,2),hn(1,1),hn(1,2),tn(1,1))
      call diff_t_testV1(n,nl,f(1,1),f(1,2),ln(1,1),ln(1,2),tn(1,2))
      call diff_t_testV1(nh,nl,hn(1,1),hn(1,2),ln(1,1),ln(1,2),tn(1,3))
      return
      end
!c-----*----------------------------------------------------6---------7--
!c     Student's t-test for the differen!ce between two ve!ctor means
!c     input: n,m,x(n,2),y(m,2),nt
!c       n: number of x
!c       m: number of y
!c       x(n,2): ve!ctor series
!c       y(m,2): ve!ctor series
!c       nt: signifi!cant level, nt must be the one of 90,95,98,99 and 999
!c     output: tn(5)
!c       tn: tn only equals 1. or 0. !corresponding to signifi!cant differen!ce or not.
!c           tn(1)~tn(5) are !corresponding to the 90%,95%,98%,99% and 99.9% !confident levels.
!c     By Dr. Jianping Li, Mar!ch 05, 2004
      subroutine diff_t_testV1(n,m,x1,x2,y1,y2,tn)
      parameter(nn=10000)
      real x(n,2),y(m,2),tn(5),x1(n),x2(n),y1(m),y2(m)
      real ft(nn,5),ax(2),ay(2),wx(n,2),wy(m,2),wx1(n),wy1(m) !work array
      do i=1,n
        x(i,1)=x1(i)
        x(i,2)=x2(i)
      enddo
      do i=1,m
        y(i,1)=y1(i)
        y(i,2)=y2(i)
      enddo
      do i=1,5
        tn(i)=0.
      enddo
      do j=1,2
        call meanvar1(n,x(1,j),ax(j))
        call meanvar1(m,y(1,j),ay(j))
      enddo
      do j=1,2
        do i=1,n
          wx(i,j)=x(i,j)-ax(j)
        enddo
        do i=1,m
          wy(i,j)=y(i,j)-ay(j)
        enddo
      enddo
      do i=1,n
        wx1(i)=sqrt(wx(i,1)**2+wx(i,2)**2)
      enddo
      do i=1,m
        wy1(i)=sqrt(wy(i,1)**2+wy(i,2)**2)
      enddo
      call meanvar(n,wx1,ax1,sx1,vx1)
      call meanvar(m,wy1,ay1,sy1,vy1)
      sxy=(n*vx1+m*vy1)/(n+m-2.)
      sxy=sxy*(1./n+1./m)
      sxy=sqrt(sxy)
      if(sxy.eq.0.)return
      d1=ax(1)-ay(1)
      d2=ax(2)-ay(2)
      dxy=sqrt(d1**2+d2**2)
      ta=dxy/sxy
      nm=n+m-2
      call t_table(ft(1,1),ft(1,2),ft(1,3),ft(1,4),ft(1,5))
      do i=1,5
        if(ta.ge.ft(nm,i))then
          tn(i)=1.
        endif
      enddo
      return
      end
!c-----*----------------------------------------------------6---------7--
!c     !computing the mean ax, standard deviation sx
!c       and varian!ce vx of a series x(i) (i=1,...,n).
!c     input: n and x(n)
!c       n: number of raw series
!c       x(n): raw series
!c     output: ax, sx and vx
!c       ax: the mean value of x(n)
!c       sx: the standard deviation of x(n)
!c       vx: the varian!ce of x(n)
!c     By Dr. LI Jianping, May 6, 1998.
      subroutine meanvar(n,x,ax,sx,vx)
      dimension x(n)
      ax=0.
      vx=0.
      sx=0.
      do i=1,n
        ax=ax+x(i)
          enddo
      ax=ax/float(n)
      do i=1,n
        vx=vx+(x(i)-ax)**2
          enddo
      vx=vx/float(n)
      sx=sqrt(vx)
      return
      end
!c-----*----------------------------------------------------6---------7--
!c     !computing the mean ax of a series x(i) (i=1,...,n).
!c     input: n and x(n)
!c       n: number of raw series
!c       x(n): raw series
!c     output: ax and sx
!c       ax: the mean value of x(n)
!c     By Dr. LI Jianping, May 6, 1999.
      subroutine meanvar1(n,x,ax)
      dimension x(n)
      ax=0.
      vx=0.
      sx=0.
      do i=1,n
        ax=ax+x(i)
      enddo
      ax=ax/float(n)
      return
      end
!c-----*----------------------------------------------------6---------7--
!c     t-distribution, i.e., student's distribution
!c     t table with two-tailed (right and left tails) probabilities
!c       P(|t|>=ta)=a
!c       where a (alpha) signifi!can!ce level and (1-a)*100% is !confinden!ce level.
!c     t90: t-distribution test at 90% !confiden!ce level.
!c     t95: t-distribution test at 95% !confiden!ce level.
!c     t98: t-distribution test at 98% !confiden!ce level.
!c     t99: t-distribution test at 99% !confiden!ce level.
!c     t999: t-distribution test at 99.9% !confiden!ce level.
!c     By Dr. Jianping Li, January 5, 2000.
      subroutine t_table(ft90,ft95,ft98,ft99,ft999)
      parameter(n=10000)
      dimension ft90(n),ft95(n),ft98(n),ft99(n),ft999(n)
      dimension n90(30),n95(30),n98(30),n99(30),n999(30)
      data n90/ 6314,2920,2353,2132,2015,1943,1895,1860,1833,1812,&
     &          1796,1782,1771,1761,1753,1746,1740,1734,1729,1725,&
     &          1721,1717,1714,1711,1708,1706,1703,1701,1699,1697/
      data n95/12706,4303,3182,2776,2571,2447,2365,2306,2262,2228,&
     &          2201,2179,2160,2145,2131,2120,2110,2101,2093,2086,&
     &          2080,2074,2069,2064,2060,2056,2052,2048,2045,2042/
      data n98/31821,6965,4541,3747,3365,3143,2998,2896,2821,2764,&
     &          2718,2681,2650,2624,2602,2583,2567,2552,2539,2528,&
     &          2518,2508,2500,2492,2485,2479,2473,2467,2462,2457/
      data n99/63657,9925,5841,4604,4032,3707,3499,3355,3250,3169, &
     &          3106,3055,3012,2977,2947,2921,2898,2878,2861,2845, &
     &          2831,2819,2807,2797,2787,2779,2771,2763,2756,2750/
      data n999/636619,31598,12941,8610,6859,5959,5405,5041,4781, &
     &     4587,4437,4318,4221,4140,4073,4015,3965,3922,3883,3850, &
     &          3819,3792,3767,3745,3725,3707,3690,3674,3659,3646/
      do i=1,30
        ft90(i)=n90(i)/1000.
        ft95(i)=n95(i)/1000.
        ft98(i)=n98(i)/1000.
        ft99(i)=n99(i)/1000.
        ft999(i)=n999(i)/1000.
      enddo
      do i=31,40
        fi=(i-30.)/10.
        ft90(i)=1.697-(1.697-1.684)*fi
        ft95(i)=2.042-(2.042-2.021)*fi
        ft98(i)=2.457-(2.457-2.423)*fi
        ft99(i)=2.750-(2.750-2.704)*fi
        ft999(i)=3.646-(3.646-3.551)*fi
          enddo
      do i=41,60
        fi=(i-40.)/20.
        ft90(i)=1.684-(1.684-1.671)*fi
        ft95(i)=2.021-(2.021-2.000)*fi
        ft98(i)=2.423-(2.423-2.390)*fi
        ft99(i)=2.704-(2.704-2.660)*fi
        ft999(i)=3.551-(3.551-3.460)*fi
          enddo
      do i=61,120
        fi=(i-60.)/60.
        ft90(i)=1.671-(1.671-1.658)*fi
        ft95(i)=2.000-(2.000-1.980)*fi
        ft98(i)=2.390-(2.390-2.358)*fi
        ft99(i)=2.660-(2.660-2.617)*fi
        ft999(i)=3.460-(3.460-3.373)*fi
  enddo
      do i=121,n
        fi=(i-120.)/(n-120.)
        ft90(i)=1.658-(1.658-1.645)*fi
        ft95(i)=1.980-(1.980-1.960)*fi
        ft98(i)=2.358-(2.358-2.326)*fi
        ft99(i)=2.617-(2.617-2.576)*fi
        ft999(i)=3.373-(3.373-2.291)*fi
  enddo
      return
      end

非我犹存

我自己有编一个标量的位势高度场的，但是不会改成矢量的。标量的如下：

        program  hecheng
    PARAMETER (Nx=144)
        PARAMETER (Ny=73)
        PARAMETER (Nt=66)
        REAL f(Nx,Ny,Nt),x(Nt),fh(Nx,Ny),fl(Nx,Ny),dh(nx,ny),dl(nx,ny),dhl(nx,ny)
        integer i,j,t,n,k,m
        real coefl,coefh
        dimension tn(nx,ny,5,3)


        OPEN(unit=15,FILE='E:\data\sdf\hgtwin500.dat',form='binary')
                read(15) (((f(i,j,t),i=1,Nx),j=1,Ny),t=1,Nt)
        close(15)

        OPEN(25,FILE='E:\data\sdf\winsn.txt')
        read(25,*) (x(t),t=1,Nt)
        close(25)

        n=66
        coefh=1
        coefl=-1
        do i=1,144
                do j=1,73                                                                          
                call diffhl1(n,x,f(i,j,:),coefh,coefl,fh(i,j),fl(i,j),dh(i,j),dl(i,j),dhl(i,j),tn(i,j,:,:))
                enddo
        enddo

        print*,'***'
        open(30,file='E:\ljp\sn\500\hgtfhsn.dat',form='binary')
        write(30) ((fh(i,j),i=1,144),j=1,73)
        close(30)
       

        open(32,file='E:\ljp\sn\500\hgtdhsn.dat',form='binary')
        write(32) ((dh(i,j),i=1,144),j=1,73)
        close(32)

        open(35,file='E:\ljp\sn\500\hgtflsn.dat',form='binary')
        write(35) ((fl(i,j),i=1,144),j=1,73)
        close(35)

        open(38,file='E:\ljp\sn\500\hgtdlsn.dat',form='binary')
        write(38) ((dl(i,j),i=1,144),j=1,73)
        close(38)
       
        open(45,file='E:\ljp\sn\500\hgtdhlsn.dat',form='binary')
        write(45) ((dhl(i,j),i=1,144),j=1,73)
        close(45)

        open(60,file='E:\ljp\sn\500\hgttnsn.dat',form='binary')
       
        do m=1,3
        do k=1,5
        do j=1,73
        do i=1,144
        write(60) tn(i,j,k,m)
        enddo
        enddo
        enddo
        enddo

        end