用nc提取出的grd和直接用nc画的图完全不同

nuistzhou

需要提取的是全球冬季感，潜热资料
利用提取的资料计算54年气候态
nc资料是高斯网格。
192*94
附件一：  用grads提取的gs

1. 'reinit'
   
2. 'sdfopen e:\bs\lhtfl.sfc.mon.mean.nc'
   
3. 'set fwrite e:\bs\global\lt.grd'
   
4. 'set gxout fwrite'
   
5. i=120  
   
6. while(i<=756)
   
7. 'set t 'i''
   
8. 'set lon 0 360'
   
9. 'set lat -88.542 88.542'
   
10. 'define lt=ave(lhtfl,t='i',t='i+2')'
    
11. 'd lt'
    
12. i=i+12
    
13. endwhile
    
14. 'disable fwrite'
    
15. 'reinit'
    
16. 'sdfopen e:\bs\shtfl.sfc.mon.mean.nc'
    
17. 'set fwrite e:\bs\global\st.grd'
    
18. 'set gxout fwrite'
    
19. i=120
    
20. while(i<=756)
    
21. 'set lon 0 360'
    
22. 'set lat -88.542 88.542'
    
23. 'define st=ave(shtfl,t='i',t='i+2')'
    
24. 'd st'
    
25. i=i+12
    
26. endwhile
    
27. 'disable fwrite'
    
28. 'reinit'
    
29. ;
    

'reinit'<br /> 'sdfopen e:\bs\lhtfl.sfc.mon.mean.nc'<br /> 'set fwrite e:\bs\global\lt.grd'<br /> 'set gxout fwrite'<br /> i=120  <br /> while(i<=756)<br /> 'set t 'i''<br /> 'set lon 0 360'<br /> 'set lat -88.542 88.542' <br /> 'define lt=ave(lhtfl,t='i',t='i+2')'<br /> 'd lt'<br /> i=i+12<br /> endwhile<br /> 'disable fwrite'<br /> 'reinit'<br /> 'sdfopen e:\bs\shtfl.sfc.mon.mean.nc'<br /> 'set fwrite e:\bs\global\st.grd'<br /> 'set gxout fwrite'<br /> i=120<br /> while(i<=756)<br /> 'set lon 0 360'<br /> 'set lat -88.542 88.542' <br /> 'define st=ave(shtfl,t='i',t='i+2')'<br /> 'd st'<br /> i=i+12<br /> endwhile<br /> 'disable fwrite'<br /> 'reinit'<br /> ;<br />

附件二：
计算气候态的f90

1. program cs
   
2. implicit none
   
3. integer,parameter::nx=192,ny=94,nt=54
   
4. real st(nx,ny,nt),lt(nx,ny,nt),rt(nx,ny,nt),ms(nx,ny),ml(nx,ny),mr(nx,ny)
   
5. integer i,j,t
   
6. open(1,file='e:\bs\global\st.grd',form='binary')               
   
7. open(4,file='e:\bs\global\mean_s.grd',form='binary')
   
8. 
   
9. do t=1,nt !nt
   
10. do j=1,ny  !ny
    
11. do i=1,nx  !nx
    
12. read(1) st(i,j,t)
    
13. enddo
    
14. enddo
    
15. enddo
    
16. 
    
17. do j=1,ny
    
18. do i=1,nx
    
19.    do t=1,nt
    
20.      ms(i,j)=ms(i,j)+st(i,j,t)
    
21.    enddo
    
22.    ms(i,j)=ms(i,j)/nt
    
23. enddo
    
24. enddo
    
25. 
    
26. do j=1,ny
    
27. do i=1,nx
    
28. write(4) ms(i,j)
    
29. enddo
    
30. enddo
    
31. 
    
32. endprogram
    

program cs<br /> implicit none <br /> integer,parameter::nx=192,ny=94,nt=54<br /> real st(nx,ny,nt),lt(nx,ny,nt),rt(nx,ny,nt),ms(nx,ny),ml(nx,ny),mr(nx,ny)<br /> integer i,j,t<br /> open(1,file='e:\bs\global\st.grd',form='binary')                <br /> open(4,file='e:\bs\global\mean_s.grd',form='binary')<br /> <br /> do t=1,nt !nt<br /> do j=1,ny  !ny<br /> do i=1,nx  !nx<br /> read(1) st(i,j,t)<br /> enddo<br /> enddo<br /> enddo<br /> <br /> do j=1,ny<br /> do i=1,nx<br />    do t=1,nt<br />      ms(i,j)=ms(i,j)+st(i,j,t)<br />    enddo<br />    ms(i,j)=ms(i,j)/nt<br /> enddo<br /> enddo<br /> <br /> do j=1,ny<br /> do i=1,nx<br /> write(4) ms(i,j)<br /> enddo<br /> enddo<br /> <br /> endprogram<br />

附件3：
      经过fotran处理后的气候态grd的描述文件

1. dset e:\bs\global\mean_s.grd
   
2. undef 1000000
   
3. xdef 192 linear  0    1.875
   
4. ydef   94 levels -88.542 -86.6532 -84.7532 -82.8508 -80.9474 -79.0435 -77.1393 -75.2351 -73.3307 -71.4262 -69.5217 -67.6171 -65.7125 -63.8079 -61.9033 -59.9986 -58.094 -56.1893 -54.2846 -52.3799 -50.4752 -48.5705 -46.6658 -44.7611 -42.8564 -40.9517 -39.047 -37.1422 -35.2375 -33.3328 -31.4281 -29.5234 -27.6186 -25.7139 -23.8092 -21.9044 -19.9997 -18.095 -16.1902 -14.2855 -12.3808 -10.476 -8.5713 -6.6666 -4.7618 -2.8571 -0.9524
   
5. 0.9524 2.8571 4.7618 6.6666 8.5713 10.476 12.3808 14.2855 16.1902 18.095 19.9997 21.9044 23.8092 25.7139 27.6186 29.5234 31.4281 33.3328 35.2375 37.1422 39.047 40.9517 42.8564 44.7611 46.6658 48.5705 50.4752 52.3799 54.2846 56.1893 58.094 59.9986 61.9033 63.8079 65.7125 67.6171 69.5217 71.4262 73.3307 75.2351 77.1393 79.0435 80.9474 82.8508 84.7532 86.6532 88.542
   
6. zdef  1 levels 1000
   
7. tdef  54 linear 00dec1957 1yr
   
8. vars 1
   
9. ms 0 99
   
10. endvars
    

dset e:\bs\global\mean_s.grd<br /> undef 1000000<br /> xdef 192 linear  0    1.875<br /> ydef   94 levels -88.542 -86.6532 -84.7532 -82.8508 -80.9474 -79.0435 -77.1393 -75.2351 -73.3307 -71.4262 -69.5217 -67.6171 -65.7125 -63.8079 -61.9033 -59.9986 -58.094 -56.1893 -54.2846 -52.3799 -50.4752 -48.5705 -46.6658 -44.7611 -42.8564 -40.9517 -39.047 -37.1422 -35.2375 -33.3328 -31.4281 -29.5234 -27.6186 -25.7139 -23.8092 -21.9044 -19.9997 -18.095 -16.1902 -14.2855 -12.3808 -10.476 -8.5713 -6.6666 -4.7618 -2.8571 -0.9524 <br /> 0.9524 2.8571 4.7618 6.6666 8.5713 10.476 12.3808 14.2855 16.1902 18.095 19.9997 21.9044 23.8092 25.7139 27.6186 29.5234 31.4281 33.3328 35.2375 37.1422 39.047 40.9517 42.8564 44.7611 46.6658 48.5705 50.4752 52.3799 54.2846 56.1893 58.094 59.9986 61.9033 63.8079 65.7125 67.6171 69.5217 71.4262 73.3307 75.2351 77.1393 79.0435 80.9474 82.8508 84.7532 86.6532 88.542<br /> zdef  1 levels 1000<br /> tdef  54 linear 00dec1957 1yr<br /> vars 1<br /> ms 0 99<br /> endvars <br />

附件4：
     用nc画的
   
此图肯定是对的

grd画的：


差异很大有没有！
求解。。。。。。。

mofangbao

楼主有没有自己先分析一下原因呢，从这么多代码中先找找原因  一行行的找，因为别人帮你找也得这么找，所以想先听听楼主自己的看法

nuistzhou

mofangbao 发表于 2012-4-23 20:17

                               
登录/注册后可看大图

楼主有没有自己先分析一下原因呢，从这么多代码中先找找原因  一行行的找，因为别人帮你找也得这么找，所以 ...

这些程序我都检查过好多遍了。没有发现问题。
我觉得最有可能是nc转grd的那个gs有问题吧

lotina

你不要作处理，使用转过来的原始grd画图，与NC画图比较，如果一致，那错就出在计算气候值的程序中，如果有错，那就是转换格式的问题了。

mofangbao

nuistzhou 发表于 2012-4-23 20:43

                               
登录/注册后可看大图

这些程序我都检查过好多遍了。没有发现问题。
我觉得最有可能是nc转grd的那个gs有问题吧

其余的没有帮你检查，至少你转换的时候不应该写 set lon 0 360  因为这两个点是同一个点，用 set x 的形式来写，否则你的图形将产生一个格点的偏移，再看你两个数据都是都是这样转换的，累计起来的当然也是错的 这个问题很早我在我空间就写过，论坛也有别人发帖子问过的。
你用lotina说的方法试一下就会知道你的数据时错的了。
另外，如果原数据中有缺测值，请在导出的grd文件中正确写上 undef，这个不能随便乱写

nuistzhou

mofangbao 发表于 2012-4-23 21:19

                               
登录/注册后可看大图

其余的没有帮你检查，至少你转换的时候不应该写 set lon 0 360  因为这两个点是同一个点，用 set x 的形式 ...

果然是set lon的问题。我是高斯网格，所以改成了 set lon 0  358.125
图基本上看起来没问题了
多谢清风！

nuistzhou

lotina 发表于 2012-4-23 21:10

                               
登录/注册后可看大图

你不要作处理，使用转过来的原始grd画图，与NC画图比较，如果一致，那错就出在计算气候值的程序中，如果有错 ...

也谢谢lotina！

lotina

nuistzhou 发表于 2012-4-23 21:47

                               
登录/注册后可看大图

果然是set lon的问题。我是高斯网格，所以改成了 set lon 0  358.125
图基本上看起来没问题了
多谢清风 ...

恭喜恭喜

hz92

很好，学习了

kxj001

学习