求海表面温度（SST）与台风生成个数的相关

小宝19950118

我用fortran求了sst与快速每年加强台风个数之间的空间相关，但是用grads画出来的图很奇怪。

下面是我的fortran程序，riyearjaso.dat存放的是每年快速加强台风的个数，jasosst.grd存放的是2×2格点上每年7-10月sst的平均值
  program laglead correlation
      parameter(nx=180,ny=89,n=66)
        real sst1(1:nx,1:ny,1:n)
        real sst(1:n)
        real x(n),y(n),rt(-1:1)
        real lag(nx,ny,-1:1)
        integer m,p,temp

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!read!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

        open(10,file='riyearjaso.dat')
        do i=1,n
        read(10,*) temp,sst(i)
        enddo



      open(20,file='jasosst.grd',form='binary')
                   do it=1,n
                do j=1,ny
                do i=1,nx
        read(20) sst1(i,j,it)
        enddo
        enddo
        enddo

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


       p=21
                do i=1,nx
                do j=1,ny
               do k=p,n
               x(k-p+1)=sst(k)
               y(k-p+1)=sst1(i,j,k)
               enddo


          nt=1
       call llcorrelation(n-p+1,x,y,nt,rt)
           do m=-1,1
          lag(i,j,m)=rt(m)
           enddo
        enddo
        enddo
       


       
          open(3,file='corr-1.grd',form='unformatted',
     $access='sequential')
             do id=-1,1
        do i=1,nx
        do j=1,ny
         write(3) lag(i,j,id)
      enddo
        ENDDO
        ENDDO
        end



c-----*----------------------------------------------------6---------7--
c     Calculating lagged and leading correlation coefficients 计算超前，滞后相关系数
c       rt(-nt:nt) between two anomaly series x(i) and y(i).
c     input: n,x(n),y(n),and nt
c       n: number of data
c       x(n) and y(n): raw anomaly series
c       nt: input parameter, maximum lag or lead time
c     output: rt(-nt:nt)
c       rt(j): lagged and leading correlation coefficients between x and y
c              j<0, lagged correlations, x trails y;
c              j=0, simultaneous correltion bewteen x and y;
c              j>0, leaing correlations, x lead y.
c     By Dr. LI Jianping, April 3, 2001.
      subroutine llcorrelation(n,x,y,nt,rt)
      dimension x(n),y(n),rt(-nt:nt)
      dimension x1(n),y1(n)               ! Work Array
c     Lagged correlations: x(i) trails y(i).
      do 20 ll=1,nt
        do 10 i=1,n-ll
          x1(i)=x(i+ll)
          y1(i)=y(i)
  10    continue
        call correlation(n-ll,x1,y1,r)
        rt(-ll)=r
  20  continue
      call correlation(n,x,y,r)
      rt(0)=r
c     Leading correlations: x(i) leads y(i).
      do 40 ll=1,nt
        do 30 i=1,n-ll
          x1(i)=x(i)
          y1(i)=y(i+ll)
  30    continue
        call correlation(n-ll,x1,y1,r)
        rt(ll)=r
  40  continue
      return
      end
c-----*----------------------------------------------------6---------7--
c     For the correlation coefficient r between two series 计算相关系数
c       x(i) and y(i), where i=1,...,n.
c     input: n,x(n),y(n)
c       n: number of time series
c       x(n): raw series
c       y(n): raw series
c     output: r
c       r: correlation coefficient between x and y
c     By Dr. LI Jianping, January 5, 2000.
      subroutine correlation(n,x,y,r)
      dimension x(n),y(n)
      call meanvar(n,x,ax,sx,vx)
      call meanvar(n,y,ay,sy,vy)
      sxy=0.
      do 10 i=1,n
        sxy=sxy+(x(i)-ax)*(y(i)-ay)
  10  continue
      sxy=sxy/float(n)
      r=sxy/(sx*sy)
      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, 1999.
      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



生成的corr-1.grd存放的是相关系数的空间分布图，用grads画图的ctl和gs文件是：
ctl
dset e:\GRADS\CORR-1\corr-1.grd
TITLE  daily  1° by 1° temperature
UNDEF  -9.99e+8
OPTIONS sequential
xdef 180 linear 0 2
ydef 89 linear  -88 2
zdef 1 levels 1000
tdef 1 linear 01JUL1970 1YR  
vars 1
corr 1 0 grid data
endvars


gs
'reinit'
'open e:\grads\corr-1\corr.ctl'
*'set x 64 68'
*'set y 70 88'
'set cmax 1'
'set cmin -1'
'set t 1'
'd corr'
;



我已经卡在这里两天了，快崩溃了。。。希望有大神可以稍微提点以下，告诉我错哪了。。。

radon

如果只是要求相关，grads有自带的函数tcorr(expr1, expr2, tdim1, tdim2)相关的说明在grads说明手册里有

小宝19950118

radon 发表于 2016-12-1 16:35
如果只是要求相关，grads有自带的函数tcorr(expr1, expr2, tdim1, tdim2)相关的说明在grads说明手册里有

嗯嗯，我试
试

小宝19950118

radon 发表于 2016-12-1 16:35
如果只是要求相关，grads有自带的函数tcorr(expr1, expr2, tdim1, tdim2)相关的说明在grads说明手册里有

谢谢！做出来了

Y.L.Zhang

路过 学习中 给楼上顶一个