求相关系数分布问题，但出不来图（又来麻烦大家了~~）

紫儿

用grads求相关系数的分布，根据论坛的提示，编辑gs文件如下，可是出不来图，并且提示如下图片，这两个ctl（下面有附）文件高度场都只有一层，但一个为0，另一个为800，所以问题如果是这里，该怎么改正才能求他们的相关系数分布
'reinit'
'open d:\abc\total.ctl'
'open d:\reason\q.ctl'
'set grads off'
'set grid off'
'set mpdset cnworld'
'set lon 100 112'
'set lat 31 43'
'set t 1 420'
*'set z 1'
'define tccn=ave(tcc,t=1,t=420)'
'define qn=ave(q.2,t=1,t=420)'
'define xy=sumg((tcc-tccn)*(q.2-qn),t=1,t=420)'
'define x2=sumg((tcc-tccn)*(tcc-tccn),t=1,t=420)'
'define y2=sumg((q.2-qn)*(q.2-qn),t=1,t=420)'
'define rr=xy/pow(x2,0.5)/pow(y2,0.5)'
'set t 1'
'set gxout shaded'
'd rr'
;

两个ctl文件如下：ctl是由nc资料编写得来，且经过验证无错，经纬一致，时间一致，分辨率一致，高度场只有一层但是高度场不一致，一个为0.一个为800
dset d:\result\netcdf.grd
title totalcloud
undef -9.99e+33
xdef 13 linear 100 1
ydef 13 linear 31 1
zdef 1 linear 0 1
tdef 420 linear 00Z01JAN1979 1mo
vars 1
tcc 0 0 Total cloud cover
endvars

dset d:\result\q.grd
title humidity
undef -9.99e+33
xdef 13 linear 100 1
ydef 13 linear 31 1
zdef 1 levels 800
tdef 420 linear 00Z01JAN1979 1mo
vars 1
q=>q 1 t,z,y,x  specific humidity
endvars

紫儿

人生路 发表于 2015-5-2 15:05
fortran怎样处理的，求大神告知

用fortran打开两个grd文件，然后按照相关系数编程
program main
parameter(nx=13,ny=13,nt=420)
integer i,j,k
real cloud(nx,ny,nt),q(nx,ny,nt),cor(nx,ny),xy(nx,ny),xx(nx,ny),yy(nx,ny),ave1(nx,ny),ave2(nx,ny)
!读入总云量要素场
open(100,file='d:/result/netcdf.grd',form='binary')
read(100)(((cloud(i,j,k),i=1,nx),j=1,ny),k=1,nt)
close(100)
!读入比湿
open(101,file='d:/result/q.grd',form='binary')
read(101)(((q(i,j,k),i=1,nx),j=1,ny),k=1,nt)
close(101)
ave1=0.0
do i=1,nx
    do j=1,ny
          do k=1,nt
       ave1(i,j)=ave1(i,j)+q(i,j,k)
          end do
          ave1(i,j)=ave1(i,j)/nt
     end do
end do
  do i=1,nx
    do j=1,ny
          do k=1,nt
      ave2(i,j)=ave2(i,j)+cloud(i,j,k)
          end do
          ave2(i,j)=ave2(i,j)/nt
      end do
   end do                                       
  xy=0
  xx=0
  yy=0
   do i=1,nx
    do j=1,ny
         do  k=1,nt
          xy(i,j)=xy(i,j)+(q(i,j,k)-ave1(i,j))*(cloud(i,j,k)-ave2(i,j))
          xx(i,j)=xx(i,j)+(q(i,j,k)-ave1(i,j))**2
         yy(i,j)=yy(i,j)+(cloud(i,j,k)-ave2(i,j))**2
         end do
        end do
end do
       
   do i=1,nx
   do j=1,ny
   cor(i,j)=xy(i,j)/(sqrt(xx(i,j)*yy(i,j)))
   end do
  end do                  
open(200,file='d:/result/corrslp.grd',form='binary')
  open(300,file='d:/result/cor.txt')
   write(200) ((cor(i,j),i=1,nx),j=1,ny)
   write(300,*) ((cor(i,j),i=1,nx),j=1,ny)
close(200)
close(300)
end program

紫儿

好可怜，怎么没有人回复我

river

你把第二个ctl里的zdef改成和第一个一样再看看
或者你使用set dfile命令把两个资料隔离开，分别计算那几个量， 最后再计算相关系数

紫儿

river 发表于 2015-5-1 07:55
你把第二个ctl里的zdef改成和第一个一样再看看
或者你使用set dfile命令把两个资料隔离开，分别计算那几个 ...

今天学校断电了，所以现在才看到，那个第一个方法试过了不行，可能是第二个z不是那个高度读不出数据的缘故，第二个办法自己查询后也不太懂那个命令，没做出来，所以能不能说详细一点，谢啦

river

你单独画第二个资料可以画出正确的图吗？

紫儿

river 发表于 2015-5-2 08:47
你单独画第二个资料可以画出正确的图吗？

可以，两个资料画出的图是对的，我现在用fortran编程处理之后再用grads画图结果出来了，但还是不明白用grads直接画图为什么不行，不过还是谢谢你的耐心解答

人生路

紫儿 发表于 2015-5-2 13:47
可以，两个资料画出的图是对的，我现在用fortran编程处理之后再用grads画图结果出来了，但还是不明白用gr ...

fortran怎样处理的，求大神告知

人生路

紫儿 发表于 2015-5-2 15:09
用fortran打开两个grd文件，然后按照相关系数编程
program main
parameter(nx=13,ny=13,nt=420)

谢谢大神

紫儿

人生路 发表于 2015-5-2 15:10
谢谢大神

我是渣，不要这么叫我，希望对你有帮助