使用的站点格点化资料不同导致Grads插值的差异

le_soleil

各位大大，我在研究5月南亚高压与我国夏季降水的相关关系的时候，使用自己的程序，所画的图如下 所使用FORTRAN处理的相关关系和转化为格点资料的程序为 rain.f90 (5.01 KB, 下载次数: 7)

2013-5-18 10:32 上传

点击文件名下载附件 ，所使用站点位置文件为，所使用gs文件为，老师的程序为，老师使用的站点位置文件为， china.map (32 Bytes, 下载次数: 3)

2013-5-18 10:38 上传

点击文件名下载附件 ，所处的图为 1.gs (618 Bytes, 下载次数: 4)

2013-5-18 10:39 上传

点击文件名下载附件 ，我发现两幅图大概还是有相似的，不过我画的值偏小，平滑，老师帮忙验证的模板并没有闭合的，但是高值区偏多，请问下问题出现在哪里啊，我不晓得是不是自己弄错了诶

le_soleil

!该程序用以计算160站资料与SAH指数的相关关系并进行显著性检验


program ex1
implicit none
integer,parameter::num=160,yt=62
external sub1
integer i,t,k,nlev,nflag
real tim,lat(num),lon(num)
character stid(num)*8
integer r6(num,yt),r7(num,yt),r8(num,yt),ave(num,yt)
real x(yt)
real Uwndx(yt),corx(num)      !中心位置指数：西风零线上H最大处的经度
real Uwndy(yt),cory(num)      !中心位置指数：西风零线上H最大处的纬度
real S(yt),cors(num)          !面积指数：(10°N-90°N,30°E-140°E)
real ERP(yt),core(num)        !东伸指数：H=M位势米环流线的最东端经度值
real I1(yt),cor1(num)         !强度指数I1：高压中心的位势高度值
real I2(yt),cor2(num)         !强度指数I2：H>=M区域内所有格点H值与M位势米差的总和
real I3(yt),cor3(num)         !强度指数I3：I2/S
real rr
!打开160站6,7,8月降水资料
open(7,file='G:\by\100hpa\1_12month\INDEX_MAY\r1606.txt')
open(8,file='G:\by\100hpa\1_12month\INDEX_MAY\r1607.txt')
open(9,file='G:\by\100hpa\1_12month\INDEX_MAY\r1608.txt')
!打开1951-2012年5月各南亚高压指数       
open(10,file='G:\by\100hpa\1_12month\INDEX_MAY\I1_5.grd',form='binary')
open(20,file='G:\by\100hpa\1_12month\INDEX_MAY\I2_5.grd',form='binary')
open(30,file='G:\by\100hpa\1_12month\INDEX_MAY\I3_5.grd',form='binary')
open(40,file='G:\by\100hpa\1_12month\INDEX_MAY\s_5.grd',form='binary')
open(50,file='G:\by\100hpa\1_12month\INDEX_MAY\uwndx_5.grd',form='binary')
open(60,file='G:\by\100hpa\1_12month\INDEX_MAY\uwndy_5.grd',form='binary')
open(70,file='G:\by\100hpa\1_12month\INDEX_MAY\erp_5.grd',form='binary')
!新建1951-2012年5月各南亚高压指数t检验
open(11,file='G:\by\100hpa\1_12month\INDEX_MAY\tI1_5.grd',form='binary')
open(21,file='G:\by\100hpa\1_12month\INDEX_MAY\tI2_5.grd',form='binary')
open(31,file='G:\by\100hpa\1_12month\INDEX_MAY\tI3_5.grd',form='binary')
open(41,file='G:\by\100hpa\1_12month\INDEX_MAY\ts_5.grd',form='binary')
open(51,file='G:\by\100hpa\1_12month\INDEX_MAY\tuwndx_5.grd',form='binary')
open(61,file='G:\by\100hpa\1_12month\INDEX_MAY\tuwndy_5.grd',form='binary')
open(71,file='G:\by\100hpa\1_12month\INDEX_MAY\terp_5.grd',form='binary')


do t=1,yt
        read(7,*)(r6(i,t),i=1,num)
        read(8,*)(r7(i,t),i=1,num)
        read(9,*)(r8(i,t),i=1,num)
end do

do t=1,yt
        read(10) I1(t)
    read(20) I2(t)
    read(30) I3(t)
    read(40) S(t)
    read(50) uwndx(t)
    read(60) uwndy(t)
    read(70) erp(t)   
end do


do i=1,num
        do t=1,yt
                ave(i,t)=r6(i,t)+r7(i,t)+r8(i,t)
        end do
end do
do i=1,num
        do t=1,yt
                ave(i,t)=ave(i,t)/3
        end do
end do

do i=1,num
        x(yt)=0.0
        do t=1,yt
                x(t)=ave(i,t)
        end do
        call sub1(x,I1,rr)
        cor1(i)=rr
    call sub1(x,I2,rr)
        cor2(i)=rr
    call sub1(x,I3,rr)
        cor3(i)=rr
    call sub1(x,S,rr)
        cors(i)=rr
    call sub1(x,uwndx,rr)
        corx(i)=rr
    call sub1(x,uwndy,rr)
        cory(i)=rr
    call sub1(x,erp,rr)
        core(i)=rr
end do



open(unit=1,file='G:\by\100hpa\1_12month\INDEX_MAY\160.txt',STATUS='OLD')
do k=1,num
    read(1,*) stid(k),lat(k),lon(k)
end do

tim=0.
nlev=1
nflag=1
do i=1,num
   write(11) stid(i),lat(i),lon(i),tim,nlev,nflag,cor1(i)
   write(21) stid(i),lat(i),lon(i),tim,nlev,nflag,cor2(i)
   write(31) stid(i),lat(i),lon(i),tim,nlev,nflag,cor3(i)
   write(41) stid(i),lat(i),lon(i),tim,nlev,nflag,cors(i)
   write(51) stid(i),lat(i),lon(i),tim,nlev,nflag,corx(i)
   write(61) stid(i),lat(i),lon(i),tim,nlev,nflag,cory(i)
   write(71) stid(i),lat(i),lon(i),tim,nlev,nflag,core(i)
end do
  
nlev = 0
write(11) stid(i-1),lat(i-1),lon(i-1),tim,nlev,nflag
write(21) stid(i-1),lat(i-1),lon(i-1),tim,nlev,nflag
write(31) stid(i-1),lat(i-1),lon(i-1),tim,nlev,nflag
write(41) stid(i-1),lat(i-1),lon(i-1),tim,nlev,nflag
write(51) stid(i-1),lat(i-1),lon(i-1),tim,nlev,nflag
write(61) stid(i-1),lat(i-1),lon(i-1),tim,nlev,nflag
write(71) stid(i-1),lat(i-1),lon(i-1),tim,nlev,nflag




close(7)
close(8)
close(9)
close(10)
close(20)
close(30)
close(40)
close(50)
close(60)
close(70)
close(11)
close(21)
close(31)
close(41)
close(51)
close(61)
close(71)
close(1)
end



subroutine sub1(x,y,r)            !子程序计算相关系数、t检验
implicit none
!integer,parameter::nt=62
!real x(nt),y(nt),rr,avex,aven,c,sx,sn
!integer it
!avex=0.0
!do it=1,nt
!        avex=avex+x(it)
!end do
!avex=avex/nt
!aven=0.0
!do it=1,nt
!        aven=aven+y(it)
!end do
!aven=aven/nt
!c=0.0
!sx=0.0
!sn=0.0
!do it=1,nt
!    c=c+(x(it)-avex)*(y(it)-aven)
!        sx=sx+(x(it)-avex)**2
!          sn=sn+(y(it)-aven)**2
!end do
!rr=c/sqrt(sx*sn)
!return
!end subroutine sub1

    integer yt,i
        parameter (yt=62)
        real x(yt),y(yt)
        real ax,ay
        real vx,vy
        real c,r
      ax=0.0
         ay=0.0
        vx=0.0
        vy=0.0
          do i=1,yt
           ax=ax+x(i)
          enddo
         ax=ax/float(yt)
       
        do i=1,yt
          ay=ay+y(i)
        enddo
         ay=ay/float(yt)
       
      
         do i=1,yt
           vx=(x(i)-ax)**2+vx
         enddo
        vx=sqrt(vx/float(yt))
       
        do i=1,yt
          vy=vy+(y(i)-ay)**2
        enddo
          vy=sqrt(vy/float(yt))

        c=0.0
           do i=1,yt
           c=(x(i)-ax)*(y(i)-ay)+c
           enddo
           c=c/float(yt)
        r=c/abs(vx*vy)
       
      return
        end
!**********主程序结束***********************
刚发现我上传的附件要收费，这是我的处理的FORTRAN程序，相关和站点转化都在一起了

le_soleil

     program correlation
        parameter (n=160,m=62)
        real ave(n,m),cor(n),con1(m),con2(m)

      open (1,file='F:\ma\i1_5.grd',FORM='binary')
      read(1) (con1(i),i=1,m)
        close(1)
        open(3,file='F:\ma\ave678.dat')
        do j=1,62
        read(3,*) (ave(i,j),i=1,160)
        enddo
         do i=1,n
         do j=1,m
         con2(j)=ave(i,j)
         enddo
!        cor(i)=coor(56,con2,con1)
         call cov(con1,con2,cor(i))
        enddo
      
        open (2,file='F:\ma\cov_r678.grd',form='binary')
        write (2) (cor(i),i=1,n)
        close (2)
        end
      subroutine cov(x,y,r)
         integer b
        parameter (b=62)
        real x(b),y(b)
        real ax,ay
        real vx,vy
        real c,r
      ax=0.0
         ay=0.0
        vx=0.0
        vy=0.0
          do i=1,b
           ax=ax+x(i)
          enddo
         ax=ax/float(b)
       
        do i=1,b
          ay=ay+y(i)
        enddo
         ay=ay/float(b)
       
      
         do i=1,b
           vx=(x(i)-ax)**2+vx
         enddo
        vx=sqrt(vx/float(b))
       
        do i=1,b
          vy=vy+(y(i)-ay)**2
        enddo
          vy=sqrt(vy/float(b))

        c=0.0
           do i=1,b
           c=(x(i)-ax)*(y(i)-ay)+c
           enddo
           c=c/float(b)
        r=c/abs(vx*vy)
       
      return
        end
!**********主程序结束***********************
这是老师的求相关的，只求了一个，进行验证的

le_soleil

    parameter(m=160)
        real vec(m)

!  read data

    open(1,file='F:\ma\cov_r678.grd',form='binary')
    READ(1)(VEC(i),i=1,m)
        close(1)
       
       
!ccccc creat grads station file

!  call stntogrd(vec1,vec2,vec3,vec4)
    call stntogrd(vec)

        end

!cccc   subroutine for transfering station-data x(160) to grid-data

    subroutine stntogrd(x)
        parameter(m=160)
    real lat(m),lon(m),x(m)
    character*8 stid(m)

!ccccc read latitude and longitude from china.dat

      open(2,file='F:\ma\CHINA.DAT',form='formatted')
      do 20 k=1,m
20   read(2,*) lat(k),lon(k)
! 12   format(f5.2,f7.2)

      do 2 i=1,m
2    stid(i)=char(i)

      OPEN(9,FILE='F:\ma\c.grd',FORM='BINARY')
       TIM=0.0
       NLEV=1
       NFLAG=1

       DO 40 I=1,m
       WRITE(9) STID(I),LAT(I),LON(I),TIM,NLEV,NFLAG,x(i)
40     continue

!cccc   On end of file write last time group terminator.
       NLEV = 0
       WRITE(9) STID(I-1),LAT(I-1),LON(I-1),TIM,NLEV,NFLAG
       close(9)
       return
       end
**************************************************
这是老师的转化为站点的程序

le_soleil

各位大大，求传教授业解惑啊

river

le_soleil 发表于 2013-5-18 10:44

                               
登录/注册后可看大图

!该程序用以计算160站资料与SAH指数的相关关系并进行显著性检验

太多了，没有谁那么有空逐个儿去给你验证，自己就不能先把范围缩小吗

le_soleil

我的绘图ctl文件
dset   G:\by\100hpa\1_12month\INDEX_MAY\terp_5.grd
dtype  station
stnmap G:\by\100hpa\1_12month\INDEX_MAY\terp.map
undef    99999.9
title    Chinese rainfall correlation datas
tdef 1  linear MAY1951 1yr
vars 1
corr  0 99 **************
endvars
我的gs文件
'reinit'
'open G:\by\100hpa\1_12month\INDEX_MAY\grid.ctl'
'open G:\by\100hpa\1_12month\INDEX_MAY\terp_5.ctl'
'set lon 70 140'
'set lat 15 55'
'set grads off'
'set grid off'
'set mpdset cnriver'
'define  cgrid=oacres(g,corr.2,10,7,4,2,1)'
'define  cgrid1=maskout(cgrid,g-0.5)'
'define  cgrid0=smth9(cgrid1)'
'set cmin 0.25'
'set gxout shaded'
'set csmooth on'
'd cgrid0'
'set gxout contour'
'set ccolor 1'
'set cint 0.1'
'cnbasemap cgrid0'
'cbarn'
'set poli off'
'run southsea -p'
'writehz 1 7 (a)erp 1  5 1 0.9 0.5 0'
'printim G:\by\100hpa\1_12month\INDEX_MAY\erp_grid.png white'
;
老师的
dset f:\ma\c.grd
dtype station
stnmap f:\ma\china.map
undef -999
title  rain and temperature  of china
tdef   1 linear jun1957 1mo
vars   1
v  0  99  summer rain of china
ENDVARS
然后是
'reinit'
'OPEN F:\ma\grid.ctl'
'OPEN F:\ma\station_cov.ctl'
'set lat 15 55'
'set lon 70 140'
'set mpdset cn'

*'DEFINE a=OACRES(G,fsr.2,10,7,4,2,1)'
'DEFINE a=OACRES(G,v.2,10,9,8,7,4,2,1)'
'DEFINE a1=MASKOUT(a,g-0.5)'
'DEFINE aa=SMTH9(a1)'
'SET GRID OFF';'SET GRADS OFF'
*'set xlopts  1 4 0.19'
*'set ylopts  1 4 0.19'
*'set xaxis 70 140 -5'
*'set yaxis 15 55 -5'
'set gxout shaded'
'set cmin 0.25'
'set csmooth on'
*'d aa'
'cnbasemap aa'
'set cthick 8'
'set clopts 1 6 0.11'
'set gxout contour'
'set cint 0.1'
'set ccolor 1'
'd aa'
'run southsea -p'

'printim F:\ma\1.png white'
;

le_soleil

river 发表于 2013-5-18 10:48

                               
登录/注册后可看大图

太多了，没有谁那么有空逐个儿去给你验证，自己就不能先把范围缩小吗

求相关的部分我已经换成一样的了，主要部分我觉得是在转化为站点上面啊，但是实在是不晓得问题会在哪里。。。。

river

le_soleil 发表于 2013-5-18 11:02

                               
登录/注册后可看大图

求相关的部分我已经换成一样的了，主要部分我觉得是在转化为站点上面啊，但是实在是不晓得问题会在哪里。 ...

你有没有把你的相关系数输出到屏幕或者txt里看看你算的相关系数对不对啊？我觉得是相关系数计算的有问题。
integer r6(num,yt),r7(num,yt),r8(num,yt),ave(num,yt)我觉得这句不能定义成integer，应该同样是实型的real。改一下试试吧。太多了，看不过来

kongfeng0824

感谢分享！！！！！！！！！！！！！！