求助！站点数据转格点场后出来的图不随t产生变化

毕业两年的小白

小白参考了清风大神的呕心之作：grads&Fortran站点绘图。站点数据写成二进制文件没有问题，可以正常出站点数据图。但在转换格点场的过程中，出现了一些问题，至今未解决，特来求助：（小白叙事比较啰嗦，还请大神见谅。真心没有随便提问啊，翻资料&查论坛搞了一天了也没解决才问的，大神勿喷）

我的数据是全国标准站的数据，仅地面一层，720*361个格点，24时次，6个变量。下面是正常做出来的站点图，站点数也足够插值了：

站点数据显示正常

站点数据显示正常

关于转换格点场，清风大神在指南中是这么说的：

小白参考大神，不考虑多时次和多层次的问题，但考虑了多变量，写了如下程序：

        integer :: i,j,k,l,m        
        real :: vargrd
        vargrd=1.0
        open(14,file='F:\temp\0718-0720-grd5.dat',status='replace',form='binary')   
            do m=1,6
                    do j=1,361
                        do i=1,720
                            write(14) vargrd
                       enddo
                   enddo
               enddo
        print*,'完成格点文件转换'
        close(14)

做出图来，无论怎么设置时间t，出来的图都一样（图片上的标题是我通过draw title手动添加的）：

小白纳闷了，又接着仔细看清风的帖子，后面还有这样一段范例：

小白没有接触过所谓的文本格式的格点数据，姑且死马当活马医，把上面的数值vargrd改为三维数组vargrd(6,361,720)一试：
    integer :: i,j,k,l,m
    real,dimension(6,361,720)  :: vargrd
    vargrd=1.0
        open(14,file='F:\temp\0718-0720-grd3.dat',status='replace',form='binary')
            do m=1,6
                    do j=1,361
                        do i=1,720
                            write(14) vargrd(m,j,i)
                     enddo
                  enddo
               enddo
        print*,'完成格点文件转换'
        close(14)

出来的图跟上面的一样，也不随t发生变化。当然这里小编很清楚，因为vargrd(m,j,i)每一个值都是1.0，与直接写入实数vargrd没有啥本质区别。
那么，是不是要加入时间和层次的循环呢？小白又把循环多套了两层：
    integer :: i,j,k,l,m
    real,       dimension(24,6,1,361,720)  :: vargrd
    vargrd=1.0
        open(14,file='F:\temp\0718-0720-grd3.dat',status='replace',form='binary')   
        do l=1,24
            do m=1,6
                do k=1,1
                    do j=1,361
                        do i=1,720
                            write(14) vargrd(l,m,k,j,i)
                     enddo
                  enddo
               enddo
            enddo
        enddo
        print*,'完成格点文件转换'
        close(14)

还是同上图，所有的图都一个模样。
之前小白还对清风大侠教程中grid数组各维的下标存在疑问，干脆把vargrd(l,m,k,j,i)换成vargrd(i,j,k,m,l)，不出所料，结果同上。


好了，现象终于说完了。有以下这么几个疑问：
①为啥我的图不随时间变化？格点场的ctl如下：
dset F:\temp\0718-0720-grd.dat
title Station Data Sample
undef 999999.0000
XDEF 720 LINEAR -180 0.5
YDEF 361 LINEAR -90 0.5
ZDEF 1 LEVELS 1000
tdef 24 linear 00z18jul2016 3hr
vars 6
ps 0 99 Surface Pressure
wd 0 99 WIN_D_Avg_10mi
ws 0 99 WIN_S_Avg_10mi
tem 0 99 tempreature
rhu 0 99 rhu
pre 0 99 pre1h
endvars

我没有用gs，就不贴了，命令与清风大侠的大同小异。
②清风大侠所说的文本格式的格点数据是什么？貌似跟我这个案例无关吧？
③如果真的需要定义多维数组vargrd，那么它的下标到底该怎么写？是vargrd（X，Y，Z，Var，T）还是倒过来？小白认为既然都是1.0，应该不影响吧？
请大神解答，研究一天了，还没弄懂/大哭
最后，贴一张小白做的彩图，也算不枉费清风大侠的苦心。效果还可以，是吧？

river

毕业两年的小白 发表于 2016-11-12 21:26
status='replace' 这个命令的意思是如果这个文件存在就重写它，如果不存在就新建，这个对结果没有影响的 ...

那我就不太清楚了，你的程序对于我只能算是看懂，所以等高手吧。如果是我的话就会把六个变量分别读到六个数组里

river

你还是把完整的写成二进制文件的程序贴出来吧

毕业两年的小白

不好意思，个人觉得站点数据可以正常出图，应该不是转换二进制的问题。不过还是贴出来源代码吧，请各位大神给诊断诊断。
代码比较啰嗦，不好意思啊。
数据就不贴了，给两张图吧：
下面这张是我代码中的 f1sta：站点经纬度表

下面这张是我代码中的f2txt：数据文件，其特点是先排列站点，每个站点的所有时次都列完后再列下一个站点的。每一行依次是：站号+时间+6个变量。


下面是我的代码：

program Catch_staID

    implicit none

!    定义变量

!    ###########################################################
!    手动输入1/2
    integer,parameter :: T=24                                     !数据文件共T个时次
    integer,parameter :: Numvar=6                             !数据文件共Numvar个变量
    integer,parameter :: X=720, Y=361, Z=1              !格点文件共有X*Y*Z个格点
    integer,parameter :: NumofStation_Id_1=2170     !站点信息文件中共有NumofStation_Id_1个站点
!    ###########################################################

    integer :: i,j,k,l,m,Line,Numsta
    real     :: tim
    integer ::nlev,flag
!    下面的变量指的是各个文件的路径
    character (len=80) :: f1sta,f2txt,f3dat,f4grd,f5ctl,f6map,f7grdctl
    character (len=1)  :: keyboard

!    下面的变量指的是站点文件
    Character*8,    dimension(NumofStation_Id_1) :: Station_Id_1
    real,            dimension(NumofStation_Id_1) :: Lat,Lon

!    下面的变量指的是数据文件，数组大小可酌情增减
    Character*8,dimension(6000)     :: Station_Id_2
    integer,    dimension(6000,4)   ::Time
    real,       dimension(6000,Numvar)  :: Var
    real,       dimension(6000)     :: Lat2,Lon2

!    下面的变量指的是格点文件
    real,       dimension(X,Y,Z,Numvar,T)  :: vargrd
    vargrd=1.0

!    ###########################################################
!    手动输入2/2：输入路径\文件名.后缀名
    f1sta='F:\My_Fortran\中国地面气象台站表 ：SURF_CHN_MUL_HOR_STATION.txt'      !站点信息文件的路径
    f2txt='C:\Users\Administrator\Desktop\T.txt'        !TXT数据文件的路径
    f3dat='C:\Users\Administrator\Desktop\T.dat'        !二进制站点数据的路径
    f4grd='C:\Users\Administrator\Desktop\T-grd.dat'    !转换后的格点数据的路径
!    ###########################################################


!    先求出数据文件中 有 效 数 据 的行数
    open(12,file=f2txt,status='old')
        Line=0
        do while(.not.eof(12))
            read(12,*)
            Line=Line+1
        enddo
        Line=Line-1
        print*,'--> 数据文件中有效数据共',Line,'行'
    close(12)

!    求数据文件站点总数
    Numsta=Line/T

!    开始一连串工作^_^

    open(11,file=f1sta,status='old')
    open(12,file=f2txt,status='old')
    open(13,file=f3dat,status='replace',form='binary')
    open(14,file=f4grd,status='replace',form='binary')

!        读取站点经纬度信息11到各个数组
        read(11,*)
        do i=1,NumofStation_Id_1
            read(11,*) Station_Id_1(i),Lat(i),Lon(i)
！            print*,  Station_Id_1(i),Lat(i),Lon(i)
        enddo
        print*, '1/6：站点信息表读取完毕'

!        读取数据文件12到各个数组
        read(12,*)
        do i=1,Line
            read(12,*) Station_Id_2(i),Time(i,:),Var(i,:)
!            print*,   Station_Id_2(i),Time(i,:),Var(i,:)
        enddo
         print*,'2/6：数据文件读取完毕'

!    检查数据文件是否有误
        i=0
        do while(i<=Line-1)
            i=i+1
            do while(Station_Id_2(i)==Station_Id_2(i+1))
                i=i+1
            enddo
            if(MOD(i,T)==0)then
                cycle
            else
                print*,'!!! 数据文件有误，站号为：',Station_Id_2(i)
                stop
            endif
        enddo
       print*,'3/6：数据文件检查完毕，行数无错误'

!        从(11)Station_Id_1中查找(12)Station_Id_2对应的经纬度
        do i=1,Line
            do j=1,NumofStation_Id_1
                if ( Station_Id_2(i)==Station_Id_1(j) ) then
                    Lat2(i)=Lat(j)
                    Lon2(i)=Lon(j)
!                    write(*,*) i,Station_Id_2(i),Lat2(i),Lon2(i)
                endif
            enddo
        enddo
        print*,'4/6：完成经纬度查找与匹配'

!        给dat文件赋值
        do i=1,T
            tim=0.0
            nlev=1
            flag=1
            do j=0,Numsta-1
                write(13) Station_Id_2(j*T+i),Lat2(j*T+i),Lon2(j*T+i),tim,nlev,flag,Var(j*T+i,:)
!                print*,   j,'-->',Station_Id_2(j*T+i),Lat2(j*T+i),Lon2(j*T+i),tim,nlev,flag,Var(j*T+i,:)
            enddo
            nlev=0
            write(13) Station_Id_2((Numsta-1)*T+i),Lat2((Numsta-1)*T+i),Lon2((Numsta-1)*T+i),tim,nlev,flag
!            print*,   j,'-->',Station_Id_2((Numsta-1)*T+i),Lat2((Numsta-1)*T+i),Lon2((Numsta-1)*T+i),tim,nlev,flag
!            pause
        enddo
       print*,'5/6：成功写入二进制文件,开始转换格点文件……'
!
!        转换出格点文件*-grd.dat
        do l=1,T
            do m=1,Numvar
                do k=1,Z
                    do j=1,Y
                        do i=1,X
                            write(14) vargrd(i,j,k,m,l)
                     enddo
                  enddo
               enddo
            enddo
        enddo
        print*,'6/6：完成格点文件转换'

    close(14)
    close(13)
    close(12)
    close(11)

end program Catch_staID

毕业两年的小白

大神在家吗？

毕业两年的小白

river 发表于 2016-11-10 21:22
你还是把完整的写成二进制文件的程序贴出来吧

大神，程序贴出来了，请您给看看

毕业两年的小白

在网上找到了阿木老师的一段授课，我的程序中是每个3h一个时刻，不知道是不是tim=0.0这里出了问题，恳请大神解答

river

毕业两年的小白 发表于 2016-11-11 21:00
不好意思，个人觉得站点数据可以正常出图，应该不是转换二进制的问题。不过还是贴出来源代码吧，请各位大神 ...

不懂最后一步那个是干嘛的？
!        转换出格点文件*-grd.dat

我觉得是你理解有误，清风的意思是假如现在你的六个变量的原始资料是格点场的文本资料，那你要写成grads可以识别的二进制文件的时候每个变量的层次循环必须分开，也就是说必须先从上到下存完一个变量，然后再存另一个变量。

而你目前制作的是一个定值为零的背景格点场，只是为了插值而已，只要最简单的直接写出来一层就够了，不需要循环时间和层次

毕业两年的小白

river 发表于 2016-11-12 17:25
不懂最后一步那个是干嘛的？
!        转换出格点文件*-grd.dat

谢谢river哥，我已经把最后一步生成格点场的代码改成了最简单的两层循环：
    integer :: i,j
    real :: vargrd
    vargrd=1.0
        open(14,file='F:\temp\0718-0720-grd.grd',status='replace',form='binary')   
                    do j=1,361
                        do i=1,720
                            write(14) vargrd
                     enddo
                  enddo
        print*,'6/6：完成格点文件转换'
        close(14)


之后做出来的图有两个问题：一是仍然不随时间变化，二是只能显示第一个变量，显示其他变量会报错：low level I/O error
所以我在此基础上添加了变量数的循环：
              do j=1,6
                    do j=1,361
                        do i=1,720
                            write(14) vargrd
                     enddo
                  enddo

              enddo
这样所有的变量都可以显示了，但是每个变量还是不随时间变化。
麻烦river哥再给看看……

如果这个问题成功解决，我一定要出一个关于多时次、多层次、多变量的帖子，算是对清风老师和您的致敬。

river

毕业两年的小白 发表于 2016-11-12 20:30
谢谢river哥，我已经把最后一步生成格点场的代码改成了最简单的两层循环：
    integer :: i,j
    rea ...

status='replace'  这个是什么意思啊？我对fortran也不太熟，但是我觉得不需要这个啊，直接就是后面form='binary' 你再试试
还有 Var(j*T+i,:) 这个也不懂，这个是什么样的数组，比如说当i=1,j=0的时候，Var(1,:)在循环过程中是否相当于第一时刻各个测站的六个变量？

毕业两年的小白

river 发表于 2016-11-12 20:49
status='replace'  这个是什么意思啊？我对fortran也不太熟，但是我觉得不需要这个啊，直接就是后面form= ...

status='replace' 这个命令的意思是如果这个文件存在就重写它，如果不存在就新建，这个对结果没有影响的。

Var(j*T+i,:) 这个数组确实像你说的那样，后面那个冒号：可以把六个变量都写出来，我通过print显示过，没有问题的。
至于这个数组前面那个下标，因为我的数据是按照站号一组一组跳的，而grads要求按照时间一组一组输入，所以需要跳着写入。