ARWpost处理后的grd数据使用fortran来读取

王磊

今天无意和朋友探讨一个问题，ARWpost处理后的WRF数据生成的grd数据，怎么用fortran语言来读取。经过一晚上的读帖子和实践，大致已经搞懂。但细细想来这样做的意义不大。既然是grd数据且配有ctl文件，最好还是应该使用grads或者将grd转化为nc数据来处理。但既然做了尝试，就姑且把代码分享出来，便于学习交流，希望对他人有所帮助。

一、问题分析


   之所以数据读取有问题，主要是在于ARWpost处理后的数据是“byteswapped”的，这从ctl文件的option选项看出，其“表示wrf模式输出数据经过ARWpost处理后，二进制的数据是反序写入的”（参见帖子：http://bbs.06climate.com/forum.php?mod=viewthread&tid=45777）。因此读入之后我们需要对数据进行一个“反序”变化，方能正确表示数据。




二、技术细节--反序处理的方法
依据数据是8字节（*8）、4字节（*4），2字节（*2）可以分为3种，但方法本质相似(这3段程序由fortran coder群提供，在此表示感谢)


（1）Elemental Function htonll( i ) result( R )  ！原始数据是8字节时使用，且处理的是整形数据

    Integer(kind=8) , Intent(IN) :: i   
    Integer(kind=8) :: R
    character :: it(8)
    it = transfer( i , it )
    R  = transfer( it(size(it):1:-1) , R )
  End Function htonll
  
(2) Ellemental Function htonl( i ) result( R )  ！原始数据是4字节时使用，且处理的是整形数据
    Integer(kind=4) , Intent(IN) :: i
    Integer(kind=4) :: R
    character :: it(4)
    it = transfer( i , it )
    R  = transfer( it(size(it):1:-1) , R )
  End Function htonl
  
（3）Elemental Function htons( i ) result( R )       ！原始数据是2字节时使用，且处理的是整形数据
    Integer(kind=2) , Intent(IN) :: i
    Integer(kind=2 :: R
    character :: it(2)
    it = transfer( i , it )
    R  = transfer( it(size(it):1:-1) , R )
  End Function htons





（3）实际演练
数据说明：我所使用的是高度场数据，ctl中的描述信息为：pdef  345 345  ；zdef   23；tdef    5，也就是说我需要一个维数为35*345*23*5的数组来存储数据，并将其转化为正常值。后续的代码如下：
PROGRAM calculate
integer,parameter :: nx=345,ny=345,nz=23,nt=5,nbite=4
real  pack_height(nx,ny,nz,nt),unpck_height(nx,ny,nz,nt)      !note the dimension order
real,external :: htonl

integer i,j,k,t

!OPEN(30,file="height.dat",access="stream")     !use ifort to compile,gfortran may not work
open(30,file="height.dat",form="unformatted",access="direct",recl=nbite*nx*ny*nz)   !use gfortran to compile, ifort may mot work
do  t=1,nt
    READ(30,rec=t)    pack_height(:,:,:,t)        !just for  gfortran
enddo
close(30)


do t=1,nt
        do i=1,nz
                do j=1,ny
                        do k=1,nx
                                unpck_height(k,j,i,t)=htonl(pack_height(k,j,i,t))   !convert the data
                                if (unpck_height(k,j,i,t) .gt. 1000000.0)  unpck_height(k,j,i,t)=-100000.0
                                !the data which greater 1000000.0 than are to be seem as missing_value,
                                !then a new one -100000.0 will be it's values
                                !the reason is that if unpck_height(k,j,i,t) is so large,
                                !so  conversion error will occaur if you output the data to a file
                        end do
                end do
        enddo
enddo


open(10,file="out.txt")
do t=1,nt
        do i=1,nz
                do j=1,ny
                        write(10,"(345(1x,f10.3))")  unpck_height(:,j,i,t)    !store the data to a file so you can  check the data
                end do
        end do
end do
close(10)
end program


Function htonl( i ) result( R )                 !依据处理方法修改，用于数据转化
    real(kind=4) , Intent(IN) :: i              !数据为实型
    real(kind=4) :: R
    character :: it(4)
        it = transfer( i , it )                       !transfer 的用法，参考fortran的学习手册即可，或者百度也能找到
    R  = transfer( it(size(it):1:-1) , R )
End Function htonl





注意：如果我们使用ifort编译，那么读取数据是改为如下格式即可，
!OPEN(30,file="height.dat",access="stream")     !use ifort to compile,gfortran may not work
do  t=1,nt
    READ(30)    pack_height(:,:,:,t)                    !just for  ifort
enddo


acess=“stream”中，stream是流的意思，以为以“数据流”形式来输入/输出文件，属于fortran2003语法？这在较新一点的编译器中才被支持。细节问题，可以参考帖子：http://v.fcode.cn/video-file_io_binary.html

夜未央

谢谢楼主！

temperature

有用有用，谢谢楼主

xugk

谢谢分享，先下了学习

哈利路亚

试用了楼主的程序，确实管用，谢谢楼主！

哈利路亚

昨天看到WRF群里doctorlee也遇到类似问题，他采用的解决方法是对文件的参数做了如下设置：
读入文件：
open(88,file=' ',form="unformatted",access="direct",recl=nbite*nx*ny*nz,convert='big_endian')
输出文件：
open(66,file=' ',form="unformatted",access="direct",recl=nbite*nx*ny*nz,convert='big_endian')
这样就不需要再对数据进行“反序”变化，直接将读入的数据write到输出文件即可。
但是不同的编译器数据格式也不尽相同，doctorlee使用以上参数设置后解决了问题，但我如此操作也不行，把输出文件的参数改为convert='little_endian',成功解决问题。
关于FORTRAN中在工作站和微机上的数据转换，给大家推荐一篇博文，相信会有启发：
http://blog.sina.com.cn/s/blog_6aabcc7d0100po05.html

对于其中具体的知识并不了解，回帖是班门弄斧了，只是把看到的解决办法贴出来，希望能对遇到相同问题的人有所帮助。

王磊

哈利路亚 发表于 2016-12-26 15:16
昨天看到WRF群里doctorlee也遇到类似问题，他采用的解决方法是对文件的参数做了如下设置：
读入文件：
op ...

感谢！刚刚我查阅了ifort编译器的手册，发现convert可以作为命令行的一个编译参数。便于数据的转换。
这样在linux上使用命令行确实可以直接转化。

王磊

这是用法，供参考。

astiny

磊啊。直接把大端改成小端不就反过来了吗，，，

xiaoguoguo327

哈利路亚 发表于 2016-12-26 15:16
昨天看到WRF群里doctorlee也遇到类似问题，他采用的解决方法是对文件的参数做了如下设置：
读入文件：
op ...

请问，你的read语句是怎么写的呢？recl的数跟read语句有关联吗？