使用GRADS处理NCEP计算水汽通量，并使用FORTRAN读取GRD文件转换为MICAPS第四类数据

yanglei_nuist

最近学习GRADS处理NCEP数据，我想先用GRADS处理60个NCEP1°×1° FNL数据，计算每天每个高度的水汽通量散度，然后输出到GRD文件中，最后使用FORTRAN读取grd文件写成MICAPS第四类格式数据，但是现在我对照GRADS直接画图和MICAPS最后呈现的图有所差异，根本对不上。
我编写的GRADS的GS文件如下：主要就是计算60个FNL里每个高度上（共21个高度）的水汽通量散度mconv，然后将每天每个高度每个格点上的计算值输出到qdiv.grd文件里，所以我认为FORTRAN里对应读取grd文件时的RECL应等于=NX*NY*4。
'reinit'
'open E:\NCEP2010\fnl_201008.ctl'
*输出到qdiv.grd中'set fwrite E:\NCEP2010\writeqdiv\qdiv.grd' m=1
while(m<61)
  'set t ' m
iz=1
  while(iz<22)
    'set x 1 360'
    'set y 1 181'
    'set z ' iz''
    'define tc=(tmpprs-273.16)'
    'define td=tc-((14.55+0.114*tc)*(1-0.01*RHprs) + pow((2.5+0.007*tc)*(1-0.01*RHprs),3) + (15.9+0.37*tc)*pow((1-0.01*RHprs),14))'
    'define vapr=6.112*exp((17.67*td)/(td+243.5))'
    'define e=vapr*1.001+(lev-100)/900*0.0034'
    'define mixr=0.62137*(e/(lev-e))*1000'
    'define qx=UGRDprs*mixr'
    'define qy=VGRDprs*mixr'
    'define mconv=hdivg(qx,qy)*1e6/9.8'
    'set gxout fwrite'
    'd mconv'
    iz=iz+1
  endwhile
  m=m+1
endwhile
'disable fwrite'
'reinit'
请注意我是一天一个高度输出一遍

然后我又写了FORTRAN文件，读取上面输出的grd文件，主要问题在读取的一步！
parameter(nx=360,ny=181,nz=21,nd=60)
integer i,j,k,l,m
real qdiv(NX,NY,NZ)
CHARACTER*2 D2(ND),T2(ND)
CHARACTER*12 DAT(ND)
CHARACTER*3 CLEV(26)
INTEGER LEV(26)
INTEGER DAY2(ND),TIM2(ND),YEAR2(ND)
CHARACTER*17 F17
YEAR2(1)=1008
DAY2(1)=10
TIM2(1)=8
WRITE(DAT(1),'(I4.4,I2.2,I2.2,A4)') YEAR2(1),DAY2(1),TIM2(1),'.000'
DO I=2,ND
        YEAR2(I)=1008
        DAY2(I)=DAY2(I-1)
        IF(MOD(I-1,4).EQ.3)THEN
                DAY2(I)=DAY2(I-1)+1
        ENDIF
        IF(MOD(I-1,4).EQ.0) THEN
                TIM2(I)=8
        ELSE IF(MOD(I-1,4).EQ.1) THEN
                TIM2(I)=14
        ELSE IF(MOD(I-1,4).EQ.2) THEN
                TIM2(I)=20
        ELSE
                TIM2(I)=2
        ENDIF
        WRITE(DAT(I),'(I4.4,I2.2,I2.2,A4)') YEAR2(I),DAY2(I),TIM2(I),'.000'
        WRITE(*,*) I,YEAR2(I),DAY2(I),TIM2(I),DAT(I)
ENDDO

LEV(1)=1000
LEV(2)=975
LEV(3)=950
LEV(4)=925
LEV(5)=900
LEV(6)=850
LEV(7)=800
LEV(8)=750
LEV(9)=700
LEV(10)=650
LEV(11)=600
LEV(12)=550
LEV(13)=500
LEV(14)=450
LEV(15)=400
LEV(16)=350
LEV(17)=300
LEV(18)=250
LEV(19)=200
LEV(20)=150
LEV(21)=100
LEV(22)=70
LEV(23)=50
LEV(24)=30
LEV(25)=20
LEV(26)=10

DO I=2,NZ
        WRITE(CLEV(I),'(I3.3)') LEV(I)
        WRITE(*,*) CLEV(I)
ENDDO

open(1,file='qdiv.grd',form='binary',recl=NX*NY*4)
do id=1,nd               
        DO K=1,NZ
                read(1,REC=(Id-1)*26+K) ((qdiv(i,j,k),i=1,NX),J=1,NY)
        ENDDO
! 这里涉及到读取grd文件，我看了fortran书，我理解recl应该等于NX*NY*4，然后REC=(Id-1)*26+K ，程序虽然能够运行，但是图片不对！
        DO I=1,NX
                DO J=1,NY
                        DO K=1,NZ
                                IF(qdiv(i,j,k)>1000.OR.qdiv(i,j,k)<-1000.0)        qdiv(i,j,k)=0.0
                        ENDDO
                ENDDO
        ENDDO

        DO IK=2,8
                OPEN(ID*100+IK,FILE='E:\NCEP2010\writeqdiv\'//CLEV(IK)//'\'//DAT(id))
                        WRITE(ID*100+IK,*) 'diamond 4 20',DAT(ID)(1:2),'年',DAT(ID)(3:4),'月',DAT(ID)(5:6),'日',DAT(ID)(7:8),'时NCEP-QDIV',CLEV(IK)
                        WRITE(ID*100+IK,*) DAT(ID)(1:2),' ',DAT(ID)(3:4),' ',DAT(ID)(5:6),' ',DAT(ID)(7:8),' ','0',' ',CLEV(IK),' ','1.0  1.0  0  359  -90  90  360  181  10 -50 50 0 0'
                        DO J=1,NY
                                WRITE(ID*100+IK,111) (qdiv(i,j,Ik),i=1,NX)
                        ENDDO
                CLOSE(ID*100+IK)
        ENDDO        

        DO IK=9,NZ
                OPEN(ID*100+IK,FILE='E:\NCEP2010\writeqdiv\'//CLEV(IK)//'\'//DAT(id))
                        WRITE(ID*100+IK,*) 'diamond 4 20',DAT(ID)(1:2),'年',DAT(ID)(3:4),'月',DAT(ID)(5:6),'日',DAT(ID)(7:8),'时NCEP-QDIV',CLEV(IK)
                        WRITE(ID*100+IK,*) DAT(ID)(1:2),' ',DAT(ID)(3:4),' ',DAT(ID)(5:6),' ',DAT(ID)(7:8),' ','0',' ',CLEV(IK),' ','1.0  1.0  0  359  -90  90  360  181  5 -20 20 0 0'
                        DO J=1,NY
                                WRITE(ID*100+IK,111) (qdiv(i,j,Ik),i=1,NX)
                        ENDDO
                CLOSE(ID*100+IK)
        ENDDO               
enddo
222 FORMAT(360(1X,F8.4))
111 FORMAT(360(1X,F8.2))

end

最后想请大家帮忙解答下，估计问题就是出现在grads输出数据或者fortran读取数据的部分（RECL的设置），谢谢大家啦！

river

好长啊，而且你说的对不上可以来个图对比一下。坐等高手，因为对micaps第四类数据格式不了解，所以也不好找问题

wuwei2163

是啊 把你的图贴上去 看看

mofangbao

建议你从你说的不对的图片找原因，把原图和你的图对比，看看错误的图形有没有什么规律，一般都会有些规律的，然后就能判断大概是哪里错了。比如经纬度反了，或者出现了偏移，或者有极大极小值等等

lqouc

说一下你的fortran用的什么编译器，有些不需要在recl里面*4的。

mofangbao

记得来反馈一下  要不然这个帖子就没意义啦

yanglei_nuist

mofangbao 发表于 2013-11-20 23:12

                               
登录/注册后可看大图

记得来反馈一下  要不然这个帖子就没意义啦

之前的那个程序也找到问题了，其实很简单，GRADS在读NCEP资料时候，我们需要设置循环，
'set fwrite E:\NCEP2010\writeqdiv\qdiv.grd'
m=1
while(m<61)
  'set t ' m
  iz=1
  while(iz<22)
    'set z ' iz''
    'define mconv=hdivg(qx,qy)*1e6/9.8'
    'set gxout fwrite'
    'd mconv'
    iz=iz+1
  endwhile
  m=m+1
endwhile
'disable fwrite'
如果这样循环，相当于把所有FNL中不同层次的所有格点上的水汽通量散度写到一个grd文件中，那么对应fortran就得这么写
open(1,file='qdiv.grd',form='binary',recl=NX*NY*4)
do id=1,nd(FNL数据个数)               
        do K=1,NZ
                read(1,REC=1) ((qdiv(i,j,k),i=1,NX),J=1,NY)！注意这里的REC一定要等于1
        ENDDO
enddo

如果像http://bbs.06climate.com/forum.php?mod=viewthread&tid=18875里的那种写GS文件
'reinit'
'open E:\NCEP2010\fnl_201008.ctl'
m=1
while(m<61)
  'set fwrite E:\NCEP2010\writeqdivday\qdiv'%m%'.grd'
  'set t ' m
  iz=1
  while(iz<22)
'set gxout fwrite'
   
    'd mconv'
    iz=iz+1
  endwhile
  'disable fwrite'
  m=m+1
endwhile
'reinit'
对应的FORTRAN就要这么写
do id=1,9       
        write(c1(id),'(i1.1)') id
        write(*,*) c1(id)
        open(id*100,file='qdiv'//c1(id)//'.grd',form='binary',recl=NX*NY*NZ*4)
        read(id*100,REC=1) (((qdiv(i,j,k),i=1,NX),J=1,NY),K=1,NZ) !这里的REC是1，而不是变量，他表示一次读多长
enddo


之前的FORTRAN程序为什么不对，就是没有理解REC的含义，REC就是表示每次读多少个RECL，上一次读一个RECL长度，下一次读2个RECL长度，那肯定是错的，所以只要合理设置REC的长度就OK了。不知道我这么解释大家能不能理解！谢谢大家的帮助啦！

guizi

楼主问题 是否 解决 期待？？？？ 我最近两天也需要捣鼓这个玩意，我感觉自己遇到的问题是access 那个选项的问题，不知道grads的输出格式 是  direct 还是  sequence 格式？？？

是否解决 期待 您的回复，我现在还不知砸个搞撒

guizi

搞定啦  我还是自己回答我的问题吧。我说我的问题，我的grads 在提取数据的时候，是顺序存放的，不知道是不是grads的提取文件都是顺序存放的，然后我读取的和你一样采用直接读取，错误时显然的，

不知道你的错误时不是在这儿，不过可以一试。希望后人在处理grads数据提取的时候，能考虑到这个环节，这个小问题，太他妈折腾哥拉 哈哈

weather3699

{:5_213:}{:5_213:}{:5_213:}