FORTRAN转MICAPS数据成grd文件 用Grads画图

安静的天分

毕业论文要求将1994-2016MICAPS地面数据相关气象要素转成grads能用的grd
现在拿1994年的出来试一下 发现stnmap这一步不行了
贴出Fortran代码 和ctl  求大神指点  
其实很早就开始弄了  不过中间有一段时间准备考研复试来着   当时拷来的数据也都有好多问题 最后把1994-2016的数据转成grd后 stnmap不了 于是就拿1994年的出来看看代码或者哪儿是不是出了问题

FORTRAN程序：
parameter(nt=1460)                                 !符号常量
integer::n1,n2,n3,n4,n,num,tim=0.0,nlev=1,nflag=1                                                                                                          
real,allocatable::lon(:),lat(:),vi(:),ac(:),ws(:)   
real,allocatable::wd(:),p(:),h(:),tt(:)            !  总云量ac，    风向wd，  风速ws， 海平面气压p，
                                                                                                        ! 6小时降水量h， 能见度vi， 温度tt,  经度lon，    纬度lat
real:: temp(2)                                     !定义了一个下标从1到2共2个元素的数组
real:: tempp(3)
real:: temppp(4)
real:: tem(1)

character*8,allocatable::sta(:)                     !动态数组，储存站号，站号必须是8个字符，5个的话会stnmap出不出来                  
character*12 filename(nt)                           !用于储存nt个时次的文件名
open(11,file='J:\MICAPSDATA\newtest\test.txt')
do i=1,nt
        read(11,*) filename(i)                      !将站号储存到filename数组里
        !print*,filename(i)                          !这是输出到最后的黑色屏幕上
enddo

close(11)
!数组filename（i）已经被写入第i个文件的 文件名
!到这儿应该没问题   

open(12,file='J:\MICAPSDATA\newtest\test.grd',status='replace',form='binary')  
!储存nt个时次的能见度场，binary表示二进制form
do k=1,nt                                           !即有nt个文件
        print*,'文件数',k                           
        open(11,file=filename(k))              
        read(11,*)
        read(11,*) n1,n2,n3,n4,n                    !将该时次的站点数赋值于n 即每个文件(即时次)里的站点数（每个文件(即时次)里的站点数不同）

                allocate(lat(n))                            !n是第k个时次（第k个000文件）的站点数目
        allocate(lon(n))
                allocate(sta(n))
        allocate(vi(n))
        allocate(ac(n))
                allocate(ws(n))
                allocate(wd(n))
                allocate(p(n))
                allocate(h(n))
                allocate(tt(n))

        print*,n                                    !从这一行往下三行 是之前一直调不通的关键 之前都直接注释掉了
        !close(11)
            !open(11,file=filename(k))                   !获知n后，重新读取。
        !read(11,*)
        !read(11,*)                                  !为的是跳过两行没有用的数据
                !现在文件指针在第三行开头了
        do i=1,n                                    !只读能见度，风速暂略        此循环为该时次的站点循环  (每个时次都有很多个站点的数据)
                                read(11,*) sta(i),lon(i),lat(i),(temp(j),j=1,2),ac(i),wd(i),ws(i),p(i),(tempp(j),j=1,3),h(i),(temppp(j),j=1,4),vi(i),(tem(j),j=1,1),tt(i)  !就是按顺序读 读完了站号经纬度然后空14格就是能见度了
                                write(12)  sta(i),lat(i),lon(i),tim,nlev,nflag,ac(i),wd(i),ws(i),p(i),h(i),vi(i),tt(i)
        enddo
close(11)

nlev=0                                             !这是一个特殊标记 告诉grads该时次结束   这个被写进了下面的那个write里
                                                   !nlev表示总层次 nlev=1表示只有地面层
                                                              !flag表示有无地面数据 有则为1
                                                           !tim=0表示地面资料 其它的暂时不管


write(12) sta(n),lat(n),lon(n),tim,nlev,nflag  !为啥有这一行先不管 至少清风大神的pdf里有这个 难道就是一个时次写完了 再把最后一个站点的前六个信息再写一下？ 其中nlev=0

!一个时次的能见度场输入完毕
deallocate(lat)
deallocate(lon)
deallocate(vi)
deallocate(sta)
deallocate(ac)
deallocate(wd)
deallocate(ws)
deallocate(p)
deallocate(h)
deallocate(tt)
!释放动态数组   
enddo
close(12)                                          
end

ctl：
DSET   J:/MICAPSDATA/newtest/test.grd
DTYPE  station
*以上用来说明这个二进制数据是一个站点格式的数据，要按站点格式的数据规范来读取
STNMAP J:/MICAPSDATA/newtest/test.map
UNDEF  9999
TITLE  test DATA
tdef 1460 linear 02Z01jan1994 6hr
vars 7
ac 0 99 all cloud
wd 0 99 wind dirtion
ws 0 99 wind speed
p  0 99 p
h  0 99 jiangshuiliang
vi 0 99 vi
tt 0 99 temperature
endvars

GRADS：
ga-> ! stnmap -i j:/MICAPSDATA/newtest/test.ctl
  Name of binary data set: J:/MICAPSDATA/newtest/test.grd
  Number of times in the data set: 1460
  Number of surface variables: 7
  Number of level dependent variables: 0

Starting scan of station data binary file.
Binary data file open: J:/MICAPSDATA/newtest/test.grd
Processing time = 1
     Time = 1 has stn count = 387
Processing time = 2
    Time = 2 has stn count = 0
Processing time = 3
  Invalid station hdr found in station binary file
  Possible causes:  Invalid level count in hdr
                    Descriptor file mismatch
                    File not station data
                    Invalid relative time
    levs = 1088421888  flag = -1064094925  time = 0

为什么第二个时次就读的0个站点呢？？？？ 求教~~