多个时次站点资料

小虹

多个时次站点资料的ctl和gs贴出来给大家参考，生成格点的fortran程序：注意红色的部分
    program main
          parameter im=3342,tm=239
          real(4) lat(im),lon(im),satdata(im),hb(im)
          real Q(im,tm)
      character*8 stid(im)
          integer station(im)
          real x(im,tm)

!read satdata
        open(100,file='E:\QC1\undef\latlon-p-final.txt')
          do i=1,im
        read(100,*) lat(i),lon(i),hb(i)
         enddo
         close(100)

        open(100,file='E:\QC1\undef\synop-q-final.txt')
        do i=1,im
        read(100,'(239f13.9)') (Q(i,j),j=1,tm)
      enddo
          close(100)
!100 format(f11.7,f9.5,f8.2)

         do i=1,im
             do j=1,tm
        x(i,j)=Q(i,j)
!         if(error(i,j).eq.0)then
!         print*,error(i,j),x(i,j)
!   x(i,j)=-999.0
!         endif
        enddo
        enddo
        open(100,file='weight-q-64.txt')
          do i=1,im
        write(100,*) lat(i),lon(i),Q(i,64)
         enddo
         close(100)

        call stntogrd(x)
        end


         subroutine stntogrd(x)
         parameter im=3342,tm=239
         real  lat(im),lon(im),satdata(im),hb(im)
         real error(im,tm)
     integer station(im)
         real x(im,tm)
         character*8 stid(im)
!READ LAT AND LON

        open(100,file='E:\QC1\undef\latlon-p-final.txt')
          do i=1,im
        read(100,*) lat(i),lon(i),hb(i)
         enddo
         close(100)



!100 format(f11.7,f9.5,f8.2)



      do  i=1,im
     stid(i)=char(i)
         enddo



!  Note: recl=11*4 means: stid*8: 2*4,  rlat,rlon tim,nlev,nflag:5*4,xyz:3*4
!  so total recl=10*4


      OPEN (3,FILE='q.grd',FORM='binary')
         do j=1,tm
                TIM=0.0
       NLEV=1
       NFLAG=1
         do i=1,im
!              if(error(i,j).eq.1)then
      WRITE(3) STID(i),LAT(i),LON(i),TIM,NLEV,NFLAG,x(i,j)

  !          print*,x(i,j)
!             endif
      enddo


       NLEV = 0
      WRITE(3) STID(im),LAT(im),LON(im),TIM,NLEV,NFLAG
          enddo
      close(3)
      end


站点ctl:
dset     E:\QC1\interploate\q-all\q.grd
dtype    station
stnmap   E:\QC1\interploate\q-all\q.map
undef    -999.0
title    observation altitude
tdef 239 linear Jan2008 1mo
vars 1
p  0 99  altitude
endvars


生成插值的格点程序：
parameter(nx=91,ny=71,nt=239)
real lat(ny),lon(nx)
real s(nx,ny,nt)
open (1,file='grid.grd',form='binary')
lat(1)=0.0         
lon(1)=55.0     
do j=1,ny-1
lat(j+1)=lat(j)+1
enddo
do i=1,nx-1
lon(i+1)=lon(i)+1
enddo
  do k=1,nt
  do i=1,nx
  do j=1,ny
   s(i,j,k)=1
   enddo
   enddo
   enddo

   write(1) s

   end


格点ctl:格点和站点时次得一致
dset E:\QC1\interploate\q-all\grid.grd
undef -999.0
xdef 91 linear 55 1
ydef 71 linear 0 1
zdef 1 linear 500 1
tdef 239 linear Jan2008 1mo
vars 1
g 0 99 pressure
endvars


铛铛铛铛，gs来了：
'reinit'
'open E:\QC1\interploate\q-all\grid.ctl'
'open E:\QC1\interploate\q-all\q.ctl'
'set grads off'
'set grid off'



'set parea 0.0 11 0.5 8'
'set grads off'
'set mpdset hires cnworld'
'set map 15 1 5'
'set map 1 1'
'set xlint 10'
'set ylint 10'
'set xlopts 1 14 0.25'
'set ylopts 1 14 0.25'

'enable print E:\QC1\interploate\q-all\q-all.gmf'

i=1
while(i<=239)
'set t 'i
'define a=oacres(g.1,p.2(t='i'))'
'define a1=maskout(a,g-0)'
'define aa=smth9(a1)'

'set gxout shaded'
'D:\study\grads-i\contourbar_qin.gs'
'set clevs  0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 8 16'
'set ccols  0 21 22 23 24 25 26 27 28 29 30 31 32 33'
'd aa*1000'
'run D:\study\grads-i\cbar_interp.gs 1 1 1 9.8 4.0 '
'basemap O 0 1'

'set line 1 1 8'
'q w2xy 55 0'
x1=subwrd(result,3)
y1=subwrd(result,6)
'q w2xy 145 70'
x2=subwrd(result,3)
y2=subwrd(result,6)

'draw rec 'x1' 'y1' 'x2' 'y2''

'draw title t='i''
'print'
*'printim E:\QC1\interploate\q-all\q-5.gif x1600 y1200 white'
'c'
i=i+1
endwhile   
'reinit'                                                                  
;


程序写的很乱，请大家忽略注释行。
多时次多变量的还有待研究

river

顶楼主！很多人都问这个问题呢，虽然乱了点儿，但有个例子总是好的。

小虹

river 发表于 2012-12-22 18:43

                               
登录/注册后可看大图

顶楼主！很多人都问这个问题呢，虽然乱了点儿，但有个例子总是好的。

谢谢~

平流层的萝卜

很给力，学习中

jeffggg

楼主 可以在每行加上注释吗？  那样好理解啊 我是新手！非常感谢

风儿飘飘

看了该贴，受益非浅！

果冻

感谢楼主的分享~受益匪浅

beiji66

单个站点多时次的资料怎么处理啊？有知道的大侠吗？

youxianchan

非常感谢~！！！

Jennifer

谢谢楼主分享好资料。