请问这个ctl和gs要怎么写啊

晓小

求谱分析的ctl和gs
c-----*----------------------------------------------------6---------7--
c     Subroutine for continuous spectrum analysis of an one-dimensional series
c         x(i) (i=1,...,n).
c     Input parameters and arrays: n,m, x(n),
c       n: number of data
c       m: biggest lag time length. Generally, m is between n/3 and n/10.
c       x(n): raw series
c     Output variables:
c        ol(0:m): frequency array
c        tl(0:m): periodic array  
c        sl(0:m): power spectrum. Its sum from 0 to m is 1.
c      st95(0:m): 95% confidence upper limit of red or white noise spectrum
c      strw(0:m): spectrum density of red or white noise
c     By Dr. Li Jianping, May 12,1998.   
       integer,parameter::n=5355,m=1100,nx=23,ny=34
         integer k,l,i
         real air(nx,ny,n)
       real x(n),sl(0:m),ol(0:m),tl(0:m),st95(0:m),strw(0:m)
       open(1,file='f:\zy\s\ecof.grd',status='old',form='binary')
       do i=1,n
          do l=1,ny
           do k=1,nx
         read(1) air(k,l,i)
        enddo;enddo;enddo
       close(1)
         print*,'1'

          do k=1,nx
           do l=1,ny
            do i=1,n
          x(i)=air(k,l,i)
            enddo
          call cspectrum(n,m,x,ol,tl,sl,st95,strw)
          enddo;enddo

         open(2,file='f:\zy\s\ps.grd',status='new',form='binary')
       do i=0,m
       write(2) ol(i),tl(i),sl(i),st95(i),strw(i)
       enddo
       print*,'2'
       stop
         end
        subroutine cspectrum(n,m,x,ol,tl,sl,st95,strw)
      dimension x(n)
      dimension sl(0:m),ol(0:m),tl(0:m),st95(0:m),strw(0:m)
      dimension r(0:m)            ! Work array
      pi=3.1415926
      call autocorrelation(n,m,x,r)
      do 10 l=0,m
        sl(l)=r(0)
        do i=1,m-1
          cc=pi*i/float(m)
          sl(l)=sl(l)+r(i)*(1.+cos(cc))*cos(l*cc)
        enddo
  10  continue
      sl(0)=sl(0)/2.
      sl(m)=sl(m)/2.
      do 20 l=0,m
        sl(l)=sl(l)/m
  20  continue
      do 30 l=1,m
        ol(l)=pi*l/m
        tl(l)=2.*m/l
  30  continue
      ol(0)=0
      tl(0)=100.*tl(1)
      call cnoisetest(n,m,r(1),sl,st95,strw)
      return
      end
c-----*----------------------------------------------------6---------7--
c     This is a subroutine for calculating the lag autocorrelations
c        r(m) of a series x(i).
c     input: n,m,and x(n)
c        n: number of the data x(n)
c        x: the raw series
c        m: maximum lagged length;  m: maximum lag time
c     output: r(0:m)
c        r: lag autocorrelations from 0 to m.
c     By Dr. LI Jianping, April 3, 1998.
      subroutine autocorrelation(n,m,x,r)
      dimension x(n),r(0:m)
      call meanvar(n,x,ax,sx,vx)
      r(0)=1.
      do i=1,n
        x(i)=x(i)-ax
      enddo
      do 10 i=1,m
        sxy=0.
        do j=1,n-i
          sxy=sxy+x(j)*x(i+j)
        enddo
        sxy=sxy/float(n)
        r(i)=sxy/vx
  10  continue
      return
      end
c-----*----------------------------------------------------6---------7--
c     Computing the mean ax, standard deviation sx
c       and variance 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 variance 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 10 i=1,n
        ax=ax+x(i)
  10  continue
      ax=ax/float(n)
      do 20 i=1,n
        vx=vx+(x(i)-ax)**2
  20  continue
      vx=vx/float(n)
      sx=sqrt(vx)
      return
      end
c-----*----------------------------------------------------6---------7--
c     Noise test for continuous spectrum (95% confidence level).
c     strw: spectrum density of red or white noise
c     chi**2 distribution with degrees of freedom v=(2*n-m/2)/m is used.
c     By Dr. Li Jianping, May 12,1998.   
      subroutine cnoisetest(n,m,r1,sl,st95,strw)
      dimension sl(0:m),st95(0:m),strw(0:m)
      dimension chi95(30),chi99(30)
      call chi_table(chi95,chi98,chi99)
      pi=3.1415926
      as=0.
      do 10 l=0,m
        as=as+sl(l)
  10  continue
      as=as/float(m+1)
      r2=r1*r1
      v=(2.*n-m/2.)/float(m)
      if(v.gt.30.)then  !! invaliable from mathmatic table for v>30
        write(*,*)'Unable continue... Please reinput m !'
        stop
      endif
      iv=int(v)
      c95=chi95(iv)+(v-iv)*(chi95(iv+1)-chi95(iv))
      if(r1.gt.0.)then
c       tname='red noise test'
        do l=0,m
          r3=1+r2-2*r1*cos(l*pi/m)
          strw(l)=as*(1.-r2)/r3
          st95(l)=strw(l)*c95/v
        enddo
      else
c       tname='white noise test'
        do l=0,m
          strw(l)=as
          st95(l)=as*c95/v
        enddo
      endif
      return
      end
c-----*----------------------------------------------------6---------7--
c     chi-square distribution
c     Chi-table with right tail probabilities
c       P(X**2>=Xa**2)=a
c       where a (alpha) significance level and (1-a)*100% is confindence level.
c     chi95: chi-square distribution at 95% confidence level.
c     chi98: chi-square distribution at 98% confidence level.
c     chi99: chi-square distribution at 99% confidence level.
c     By Dr. Li Jianping, May 12,1998.   
      subroutine chi_table(chi95,chi98,chi99)
      dimension chi95(30),chi98(30),chi99(30)
      dimension c95(30),c98(30),c99(30) !work array
      data c95/ 3.841, 5.991, 7.815, 9.488,11.070,
     *         12.592,14.067,15.507,16.919,18.307,
     *         19.675,21.026,22.362,23.685,24.996,
     *         26.296,27.587,28.869,30.144,31.410,
     *         32.671,33.924,35.172,36.415,37.652,
     *         38.885,40.113,41.337,42.557,43.773/
      data c98/ 5.412, 7.824, 9.837,11.668,13.388,
     *         15.033,16.622,18.168,19.679,21.161,
     *         22.618,24.054,25.472,26.873,28.259,
     *         29.633,30.995,32.346,33.687,35.020,
     *         36.343,37.659,38.968,40.270,41.566,
     *         42.856,44.140,45.419,46.693,47.662/
      data c99/ 6.635, 9.210,11.345,13.277,15.086,
     *         16.812,18.475,20.090,21.666,23.209,
     *         24.725,26.217,27.688,29.141,30.578,
     *         32.000,33.409,34.805,36.191,37.566,
     *         38.932,40.289,41.638,42.980,44.314,
     *         45.642,46.963,48.278,49.588,50.892/
      do i=1,30
        chi95(i)=c95(i)
        chi98(i)=c98(i)
        chi99(i)=c99(i)
      enddo
      return
      end

言深深

好长··········

晓小

从李建平老师主页上下的