grads画站点图老是有问题

林宥藜

    附上程序：（附件是所需要的.ctl .gs .grd文件）
  program EP
      implicit none
      integer,parameter::nyr=61,n=160,nx=144,ny=73
real a(160,60),ind(61),r(160),aa(61),bb(61),lat(160),
     &     lon(160),h(nx,ny,12,nyr),eu(nyr),aveeu,sxeu,aeu(61),ssxeu
      real aveh(nx,ny),sxh(nx,ny),reh(nx,ny),avet(n),sxt(n),aaveu
      real tim,nflag,nlev
!eu--遥相关，aeu--标准遥相关，sxeu--方差，sxh--h的方差，reh--遥相关和高度场的相关系数
!ssxeu--58年的eu方差，aaveu--58年eu平均（1951-2008）
character*8 id(160)
      integer i,k,j,it
open(2,file='d:\dqsj03\data\t1601.txt')
open(3,file='d:\dqsj03\ind.grd',form='binary')
open(4,file='d:\dqsj03\data\lat_lon.txt')
      open(5,file='d:\dqsj03\r160.grd',form='binary')
      open(1,file='d:\dqsj03\data\hgt500.grd',form='binary')
      open(19,file='d:\dqsj03\r160.dat')
ccccccccccccccc 读数据（指数、经纬度、160站温度）
ccc    a:160站气温(1951~2010年)   ind：指数序列(1948~2008年)
      read(2,*),((a(i,j),i=1,160),j=1,60)
      do i=1,160
      do j=1,60
      a(i,j)=a(i,j)/10
      enddo;enddo
      do i=1,160
      read(4,*),lat(i),lon(i)
      enddo
      do k=1,nyr
      do it=1,12
      do j=1,ny
      do i=1,nx
      read(1) h(i,j,it,k)
      enddo; enddo; enddo;enddo
ccccccccccccccc 编程求相关
      open(7,file='d:\dqsj03\data\eu.grd',form='binary')
      open(8,file='d:\dqsj03\data\aeu.grd',form='binary')
      open(9,file='d:\dqsj03\data\reh.grd',form='binary')
      do k=1,61
      eu(k)=-0.25*h(9,59,1,k)+0.5*h(31,59,1,k)-0.25*h(59,53,1,k)
      enddo
      do k=1,61
      write(7),eu(k)
      write(3),eu(k)
      enddo
      !print *,eu(1)
      aveeu=0.0
      sxeu=0.0
      do k=1,nyr
      aveeu=aveeu+eu(k)
      enddo
      aveeu=aveeu/nyr
       !print *,aveeu
      do k=1,nyr
      sxeu=(eu(k)-aveeu)*(eu(k)-aveeu)+sxeu
      enddo
       sxeu=sxeu/61
       !print *,sxeu
      do k=1,61
      aeu(k)=(eu(k)-aveeu)/sqrt(sxeu)
      enddo
      do k=1,61
      !print *,aeu(i)
      write(8),aeu(k)
      enddo
      aveh(nx,ny)=0.0
      do j=1,ny
      do i=1,nx
      do k=1,61
      aveh(i,j)=h(i,j,1,k)+aveh(i,j)
      enddo
      aveh(i,j)=aveh(i,j)/nyr
      !print *,aveh(i,j)
      enddo;enddo
      sxh(nx,ny)=0.0
      do j=1,ny
      do i=1,nx
      do k=1,61
      sxh(i,j)=(h(i,j,1,k)-aveh(i,j))*(h(i,j,1,k)-aveh(i,j))+sxh(i,j)
      enddo
      sxh(i,j)=sxh(i,j)/61
      !print *,sxh(i,j)
      enddo;enddo
      reh(nx,ny)=0.0
      do j=1,ny
      do i=1,nx
      do k=1,61
      reh(i,j)=reh(i,j)+(eu(k)-aveeu)*(h(i,j,1,k)-aveh(i,j))
      enddo
      reh(i,j)=reh(i,j)/sqrt(sxeu*sxh(i,j))/61
      write(9),reh(i,j)
      enddo;enddo
      avet=0.0
      do i=1,160
      do j=1,58
      avet(i)= avet(i)+a(i,j)
      enddo
      avet(i)=avet(i)/58
      !print *,avet(i)
      enddo
      sxt(n)=0.0
      do i=1,n
      do j=1,58
      sxt(i)=(a(i,j)-avet(i))*(a(i,j)-avet(i))+sxt(i)
      enddo
       sxt(i)=sxt(i)/58
      !print *,sxt(I)
      enddo
      aaveu=0.0
      ssxeu=0.0
      do k=4,nyr
      aaveu=aaveu+eu(k)
      enddo
      aaveu=aaveu/58
      do k=4,nyr
      ssxeu=(eu(k)-aaveu)*(eu(k)-aaveu)+ssxeu
      enddo
      ssxeu=ssxeu/58
      !print *,ssxeu
      r(n)=0.0
      do i=1,n
      do j=1,58
      r(i)=r(i)+(eu(j+3)-aveeu)*(a(i,j)-avet(i))
      enddo
      r(i)=r(i)/sqrt(sxt(i))/sqrt(ssxeu)/58
         !print *,ret(i)
      enddo
ccccccccccccccccccc写站点数据
  do j=1,160
  id(j)=char(j)
   tim=0.0
  nlev=1
  nflag=1
  write(5)id(j),lat(j),lon(j),tim,nlev,nflag,r(j)
      write(19,*) id(j),lat(j),lon(j),tim,nlev,nflag,r(j)
      enddo
  tim=0.0
  nlev=0
  nflag=1
  write(5)id(j-1),lat(j-1),lon(j-1),tim,nlev,nflag
      write(19,*) id(j-1),lat(j-1),lon(j-1),tim,nlev,nflag
      end
映射文件一直不能生成：
[img]file:///C:/Users/ASUS/AppData/Roaming/Tencent/Users/826499352/QQ/WinTemp/RichOle/FWG4JM3BUF(Q(~[9WFPUT5I.png[/img]

林宥藜

dset d:\dqsj03\r160.grd
dtype station
stnmap d:\dqsj03\r160.map
undef 99999.9
title the 160 station Jan t and ind rc
tdef 1 linear jan1984 1yr
vars 1
r  0 99 winter rain anomaly
endvars
;

林宥藜

.gs
'reinit'
'open d:\dqsj03\data\map.ctl'
'open d:\dqsj03\r160.ctl'

'set lat 15 55'
'set lon 70 140'
'set xlopts 1 6 0.15'
'set ylopts 1 6 0.15'
'define aa=oacres(g.1,r.2)'
'define rgg=maskout(aa,g.1-0.5)'
'define bb=smth9(rgg)'
'set grads off'
'set mpdset lowres'
'set poli on'
'set csmooth off'
'set mpdraw on'
'set map 1 1 1'

'set gxout shaded'
'set clevs -0.32 -0.25 0.25 0.32'
'set ccols 2 5 0 9 3'
'd bb'
'set gxout contour'
'set cint 0.1'
'set gxout contour'
'set cthick 3'
'd bb'
'set cthick 6'
'set clevs 0'
'set ccolor 1'
'd bb'
'run d:\dqsj03\map\china.gs'
'run d:\dqsj03\map\taiwan.gs'
'run d:\dqsj03\map\hainan.gs'
'run d:\dqsj03\map\river.gs'
'run d:\dqsj03\map\southsea.gs'

'enable print d:\dqsj03\r160.gmf'
'print'
'disable print'
;