Grads多个时次循环输出数据到同一个文件，只有第一个时次数据能正常作图

银酱赛高

grads读取nc数据，计算整层水汽通量积分，多个时次循环并输出到同一个文件出现问题
1、首先，附上数据ctl：

2、然后，计算整层水汽通量，垂直层次选取1000hPa-100hPa，挑选11个时次，输出每个时次选定区域的整层水汽通量积分，为grd格式，gs文件如下：

1. 'reinit'
   
2. 'sdfopen h:\ecmwf_typhoon\200806.nc'
   
3. 
   
4. 'set lat 15 30'
   
5. 'set lon 105 125'
   
6. 'set z 1'
   
7. 'set t 101 111'
   
8. 
   
9. 'define e=exp(17.67*(t-273.15)/(t-29.65))*(r/100)*6.112'
   
10. 'define q=0.622*e/(lev-0.378*e)'
    
11. 'define qu=u*q/9.8'
    
12. 'define qv=v*q/9.8'
    
13. 'define quall=vint(z(z=1),qu,100)'
    
14. 'define qvall=vint(z(z=1),qv,100)'
    
15. 
    
16. 'set fwrite h:\ecmwf_typhoon\quvall.grd'
    
17. 'set gxout fwrite'
    
18. 'set lat 15 30'
    
19. 'set lon 105 125'   
    
20. 
    
21. t1=101
    
22. while(t1<=111)
    
23.     'set t 't1''
    
24.     'set z 1'
    
25.     'd quall'
    
26.     'd qvall'  
    
27.     t1=t1+1
    
28. endwhile
    
29.    
    
30. 'disable fwrite'
    
31. ;

复制代码

3、运行gs文件没有报错，为grd文件写ctl文件如下：

4、‘set t 1’，'d quall'能正常出图，但是t=2，……11出图都有问题。


5、排查代码问题，一开始怀疑是多时次循环数据写入文件字段有问题，尝试循环输出其他变量（气温），该段代码更改如下：

1. t1=101
   
2. while(t1<=111)
   
3.     'set t 't1''
   
4.     'set z 1'
   
5.     'd t'
   
6.     t1=t1+1
   
7. endwhile

复制代码

发现气温能够正常输出并作图，证明该段代码没有错误，ctl描述文件也没问题。
6、然后猜测可能是整层水汽通量积分计算方法有问题，尝试逐时次输出整层水汽通量的分布图，代码更改如下：

1. t1=101
   
2. while(t1<=111)
   
3.     'set t 't1''
   
4.     'set z 1'
   
5.     'set gxout contour'
   
6.     'd quall'
   
7.     'd qvall'  
   
8.     'printim h:\ecmwf_typhoon\'%t1%'.png white'
   
9.     'c'
   
10.     t1=t1+1
    
11. endwhile

复制代码

发现每个时次均可以正常出图，证明计算方法代码也没有错：
求问各位大神，gs文件怎样修改才能每个时次的数据都能输出并正常作图？？

苏拉苏拉

ctl缺测设为-9.99+08试试

river

1. 'reinit'
   
2. 'sdfopen h:\ecmwf_typhoon\200806.nc'
   
3. 
   
4. 'set gxout fwrite'
   
5. 'set fwrite h:\ecmwf_typhoon\quvall.grd'
   
6. 'set lat 15 30'
   
7. 'set lon 105 125'   
   
8. 
   
9. t1=101
   
10. while(t1<=111)
    
11.     'set t 't1''
    
12.     'set z 1'
    
13. 'define e=exp(17.67*(t-273.15)/(t-29.65))*(r/100)*6.112'
    
14. 'define q=0.622*e/(lev-0.378*e)'
    
15. 'define qu=u*q/9.8'
    
16. 'define qv=v*q/9.8'
    
17. 'define quall=vint(z(z=1),qu,100)'
    
18. 'define qvall=vint(z(z=1),qv,100)'
    
19. 'd  quall'
    
20. 'd  qvall'
    
21.     t1=t1+1
    
22. endwhile
    
23.    
    
24. 'disable fwrite'
    
25. ;
    
26. 
    
27. 
    

复制代码

银酱赛高

苏拉苏拉 发表于 2020-6-6 13:42
ctl缺测设为-9.99+08试试

可以了！台风中心的整层水汽通量积分存在极值中心，一开始忽略了这点。谢谢！

银酱赛高

river 发表于 2020-6-6 23:51

问题已经解决了，谢谢您