GrADS站点资料作图详细解决方案

风儿飘飘

将错误的信息发上来，请给看看

mofangbao

风儿飘飘 发表于 2013-8-18 20:54

                               
登录/注册后可看大图

将错误的信息发上来，请给看看

一般是生成二进制文件的时候循环上出错了，没有完全按照规范些，当然你在没找到错误之前会觉得你是按照规范些的~

风儿飘飘

!************************变量说明**********************************
! stn：站点数（固定值），year：年份，固定参数，每年修改一次即可
! lon/lat：经纬度数组
! x：读取出来的某月160数据
! stid：站号
! lab：标记要转换的是降水还是温度，kind：温度降水字符串
! imon,iyea：整型的月份和年份
! mon,yea：字符型的月份和年份（用于生成文件名）
! filename：打开的文件名称（由输入的年份和月份而定）
! bt：生成的文件名称
! next：记录是否继续转换
!资料从1951年开始
!************************************************************************************

     program ex_03
      integer,parameter::stn=160,year=62
      integer::i,j,k
      real ::lat(stn),lon(stn)
      real ::x(stn,year),y(stn,year),z(stn,year)
     real ::xyz(stn,year),xyzmean(stn,year)
      integer  ::stid(stn)
     real ::tim
    integer ::lev,flag
!!!!!!!!!!!!!read 160 station precipitation!!!!!!!!!!!!!!!!!!!
      print*,'begin to read precipitation'
      open(1,file="F:/China160station/precipitation/r1606.txt",
     $        status='old')
      do j=1,year
      read(1,*)(x(i,j),i=1,stn)
      enddo
      close(1)
        print*,'read June precipitation over'

        print*,'begin to read precipitation'
      open(1,file="F:/China160station/precipitation/r1606.txt",status='old')
      do j=1,year
      read(1,*)(y(i,j),i=1,stn)
      enddo
      close(1)
        print*,'read July precipitation over'

        print*,'begin to read precipitation'
      open(1,file="F:/China160station/precipitation/r1606.txt",status='old')
      do j=1,year
      read(1,*)(z(i,j),i=1,stn)
      enddo
      close(1)
        print*,'read August precipitation over'

        do j=1,year
        do i=1,stn
        xyz(i,j)=x(i,j)+y(i,j)+z(i,j)
        xyzmean(i,j)=xyz(i,j)/3
        enddo
        enddo



      
!!!!!!!!!!!!!!!!!zhuanhuanrejizhi!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!111
      print*,'begin to read 160 stations'
      open(2,file='F:/China160station/precipitation/station160.txt',status='old')
      do k=1,stn
      read(2,*)stid(k),lat(k),lon(k)
!        print*,stid(k),lon(k),lat(k)
!        pause
      enddo
      close(2)
        print*,'successfully write station data'

  !!!!!!!!!!!!!!!!!!按顺序写!!!!!!!!!!!!!!!!!
open(4,file='F:/China160station/sta2grd/sta.grd',status='replace', form='binary')
        do j=1,year
                 tim=0.0
        nlev=1
        flag=1
        do i=1,stn
        write(4)stid(i),lat(i),lon(i),tim,nlev,flag,xyzmean(i,j)
        enddo
        nlev=0
        write(4)stid,lat,lon,tim,nlev,flag
             enddo
        close(4)
        print*,'successfully write data'
        end

风儿飘飘

mofangbao 发表于 2013-8-18 22:33

                               
登录/注册后可看大图

一般是生成二进制文件的时候循环上出错了，没有完全按照规范些，当然你在没找到错误之前会觉得你是按照规 ...

谢谢清风，我把读取160站fortran程序贴上来了，就是最后输出到文件4时有问题，但反复看了，不知问题出在哪，我看了您的贴子，不知如何解决，请给点意见,谢谢！

mofangbao

风儿飘飘 发表于 2013-8-19 12:37

                               
登录/注册后可看大图

谢谢清风，我把读取160站fortran程序贴上来了，就是最后输出到文件4时有问题，但反复看了，不知问题出在哪 ...

write(4)stid,lat,lon,tim,nlev,flag
你这样会把数组中所有数据都输出的，最好理解了再复制粘贴

风儿飘飘

mofangbao 发表于 2013-8-19 14:26

                               
登录/注册后可看大图

write(4)stid,lat,lon,tim,nlev,flag
你这样会把数组中所有数据都输出的，最好理解了再复制粘贴

不理解就是为什么nlev=0时，表示一个时次的输出结束，为何还要再写write(4)stid,lat,lon,tim,nlev,flag，输出下一时次的不是进行下一个循环吗，请给点意见

mofangbao

风儿飘飘 发表于 2013-8-20 09:07

                               
登录/注册后可看大图

不理解就是为什么nlev=0时，表示一个时次的输出结束，为何还要再写write(4)stid,lat,lon,tim,nlev,flag， ...

以后遇到这种不理解的问题希望第一反应是去官方找参考文档哈：
http://www.iges.org/grads/gadoc/aboutstationdata.html#station

风儿飘飘

谢谢清风！

dannyjoyride

太感谢LZ了，正在学习中。。。

cx0515

非常感谢，非常感谢，很实用的