求助：学习利用micaps地面资料画6h降水图时遇到的问题

F.Lee

第一次做micaps多时次单变量资料的批处理，学习了论坛里好多大神的帖子，借鉴了luoziwuhui的帖子，http://bbs.06climate.com/forum.php?mod=viewthread&tid=13648，编写fortran程序读取数据并且转成二进制文件，但是在程序运行的过程中总是遇到如图的问题
我用的是micaps   surface文件夹下r6-p文件夹下的文件，下面是数据截图

我用Windows 的CMD命令将要用的17个时次的数据写入了filename.txt，不知道为什么写出来不是像别人一样每行一个文件名，而是每行三个，由于文件数比较少，我就手动改了。


fortran程序如下：
program read_surface_observation
          parameter(nt=17)
          integer nlev,flag,n1,n
           real,allocatable::lon(:),lat(:),rd(:),var(:)    !动态数组分别用来储存站点的降水
          real tim
                  character*8,allocatable::sta(:)      !动态数组，储存站号，站号必须是8个字符，5个的话会stnmap出不出来   
        character*12 filename(nt)        !用于储存32个时次的文件名
           open(1,file='G:\rain\filename.txt')
           do i=1,nt
               read(1,*) filename(i)
               print*,'Filename:',filename(i)
            enddo
            close(1)               
        open(2,file='G:\rain\rain6h.grd',status='replace',form='binary')  !储存32个时次的降水
       do k=1,nt
        open(3,file=filename(k))
        read(3,*)
                read(3,*)
                read(3,*)
                read(3,*)
                read(3,*)
                read(3,*)
                read(3,*)
                read(3,*)
                read(3,*)
                read(3,*)
                read(3,*)
                read(3,*)
                read(3,*)
        read(3,*) n1,n    !将该时次的站点数赋值于n
        allocate(lat(n))  !分配动态数组
        allocate(lon(n))
                allocate(var(n))
        allocate(rd(n))
        allocate(sta(n))
        close(3)
        open(4,file=filename(k))   !获知n后，重新读取数据。
        read(4,*)
        read(4,*)
          do i=1,n        !只读降水，别的要素暂略
             read(4,*) sta(i),lon(i),lat(i),var(i),rd(i)
              tim=0.0;nlev=1;nflag=1
              write(2)  sta(i),lat(i),lon(i),tim,nlev,flag,rd(i)  !一个时次的降水场输入完毕
           enddo
           nlev=0
            write(2) sta(n),lat(n),lon(n),tim,nlev,flag
        !一个时次的降水场输入完毕，告诉Grads 该时次的数据结束。这是一个特殊标记。
        deallocate(lat)
        deallocate(lon)
        deallocate(rd)
        deallocate(sta)
!释放动态数组
        close(4)
        enddo
         close(2)
      end program read_surface_observation

麻烦学姐学长帮忙看看，哪里出了问题，我已经尽量把问题描述清楚了，可能你们会觉得我的问题过于小白，借鉴人家的都弄不懂，可是我确实很努力了，纠结这个问题好久了，身边的同学也没有会的。我在原帖下面问楼主，可能楼主比较忙，没理我。。。真心希望有人能帮助我一下，可能你的一句话就会让我柳暗花明。学妹在这里先谢过了！

平流层的萝卜

这个程序我当时写的很臃肿，没必要用动态数组的。一堆read（3，*）也可以用循环代替。回头帮你看一下。先吃饭了。

F.Lee

平流层的萝卜 发表于 2014-4-11 18:08
这个程序我当时写的很臃肿，没必要用动态数组的。一堆read（3，*）也可以用循环代替。回头帮你看一下。先吃 ...

先多谢了！！！

lightxxia

可能是变量定义的问题 你一开始定义的站号是字符型的 可是你读取的站点数是整形的 在写入之前把站点数改成字符型的试下把 用sta（i）=char（st（i）） 我之前也遇到类似的问题 最后可以了 好像我们处理的是同样的数据 你合成的时候对了么 可以的话我们可以交流下 我现在只是单独的分析了每个时次的 还没有合成后分析的

平流层的萝卜

先说说简化一下程序吧，因为当时写的代码实在臃肿，很难看。1,首先删除有关动态数组的声明及调用。因为sta(i),lat(i),lon(i),tim,nlev,flag,rd(i)完全没必要用数组，在获取站数的前提下，用站数循环，循环体内直接用变量sta,lat,lon,tim,nlev,flag,rd就好，省去了不少麻烦事。

2,13个连续的    read(3,*)用
do i=1,13
    read(3,*)
enddo

代替
3.无需open(3,file=filename(k))。。。。 。close(3) 再open(4,file=filename(k)) ，直接一次打开该文件，在里边先获取n站数，然后读取站点降水，写入grd。
再说你的错误。
open(4,file=filename(k))   !获知n后，重新读取数据。
read(4,*)
read(4,*)
后边接的是读取站数，为什么要这么写呢？你的站点将水是从第15行开始的。你直接读取14行的站数，再依次向下读取站点将水值就好了啊。你的意思是从第3行开始读取的。不理解。

平流层的萝卜

lightxxia 发表于 2014-4-11 21:50
可能是变量定义的问题 你一开始定义的站号是字符型的 可是你读取的站点数是整形的 在写入之前把站点数改成 ...

站点数n是整型的没错啊，你表达的意思是站点名称改成字符型的吧，程序里已经把站点名称定义成字符型的了。这样读取是可以的。

F.Lee

lightxxia 发表于 2014-4-11 21:50
可能是变量定义的问题 你一开始定义的站号是字符型的 可是你读取的站点数是整形的 在写入之前把站点数改成 ...

插入一句sta(i)=char(st(i))之后出现错误和警告了，
G:\rain\for\plot.f90(43) : Warning: Arguments' data types are incompatible with intrinsic procedure, assume EXTERNAL.   [CHAR]
     sta(i)=char(st(i))
---------------------------------^
G:\rain\for\plot.f90(43) : Error: A CHARACTER data type is required in this context.   [CHAR]
     sta(i)=char(st(i))

F.Lee

平流层的萝卜 发表于 2014-4-11 22:15
先说说简化一下程序吧，因为当时写的代码实在臃肿，很难看。1,首先删除有关动态数组的声明及调用。因为sta( ...

感谢前辈！read(4,*)确实写错了，可是我改成了从第15行开始读取降水值，问题依然存在，怎么办呢？

平流层的萝卜

再提供一种可能的解决方案，你先试试，在你的1520.000文件最后一行的末尾，敲击回车。再试试。

F.Lee

平流层的萝卜 发表于 2014-4-11 22:33
再提供一种可能的解决方案，你先试试，在你的1520.000文件最后一行的末尾，敲击回车。再试试。

还是不行啊