MATLAB计算假相当位温的程序，请大神帮忙检查一下有没有错误。谢谢！

故乡月(环)

利用坛子里面计算假相当位温的GS文件，自己修改了一下，绘制东经103-109°，北纬26-30°区域的平均假相当位温剖面图，不知有错否，请大神指点。资料为1°*1°的fnl资料，直接下载netcdf格式。

1. 
   
2. % function [ output_args ] = eqtdraw( input_args )
   
3. %EQTDRAW Summary of this function goes here
   
4. %   Detailed explanation goes here
   
5. %   计算假相当位温并画图
   

复制代码

故乡月(环)

1. 
   
2. % function [ output_args ] = eqtdraw( input_args )
   
3. %EQTDRAW Summary of this function goes here
   
4. %   Detailed explanation goes here
   
5. %   计算假相当位温并画图
   
6. 
   
7. clc;clear all
   
8. ncdisp('D:\fnl\fnl_20150418_12_00.grib2.nc');
   
9. prs=(ncread('D:\fnl\fnl_20150418_00_00.grib2.nc','lv_ISBL0'))/100
   
10. t1=ncread('D:\fnl\fnl_20150418_00_00.grib2.nc','TMP_P0_L100_GLL0');
    
11. rh1=ncread('D:\fnl\fnl_20150418_00_00.grib2.nc','RH_P0_L100_GLL0');
    
12. 
    
13. 
    
14. lon=ncread('D:\fnl\fnl_20150418_12_00.grib2.nc','lon_0');
    
15. lat=ncread('D:\fnl\fnl_20150418_12_00.grib2.nc','lat_0');
    
16. S2=shaperead('D:\国家基础地理信息系统数据\国界与省界\bou2_4p.shp');
    
17. bou1x=[S2(:).X];
    
18. bou1y=[S2(:).Y];
    
19. 
    
20. for j=1:26
    
21.    lev(1:360,1:181,j)=prs(j);
    
22. end
    
23. 
    
24. es1=(6.112.*exp(17.67.*(t1-273.15)./(t1-29.65)));
    
25. q1=rh1.*(0.62197.*es1./(lev-es1))/100;
    
26. e1=lev.*q1./(0.62197+q1)+1e-10;
    
27. tlcl1=55.0+2840./(3.5.*log(t1)-log(e1)-4.805);
    
28. theta1=t1.*(1000./lev).^(0.2854.*(1.0-0.28.*q1));
    
29. eqt1=theta1.*exp(((3376./tlcl1)-2.54).*q1.*(1.0+0.81.*q1))-273.5;
    
30. tt1=eqt1(103:109,61:65,:);
    
31. gg1=mean(mean(tt1));
    
32. gp1=squeeze(squeeze(gg1));
    
33. plot(flipud(gp1(6:end)),flipud(prs(6:end)),'-k*')
    
34. set(gca,'YDir','reverse')
    
35. axis([50 120 100 1000])
    
36. % print('-dpng','D:\7.png','-r300')

复制代码

yuppieflu

看了一下基本没错的，楼主所用公式是否摘自Bolton，1980？还有，e1=lev.*q1./(0.62197+q1)+1e-10;这一条+1e-10是什么意思？

呼啦啦

e1=lev.*q1./(0.62197+q1)+1e-10；中最后的1e-10是什么意思？

呼啦啦

yuppieflu 发表于 2016-12-5 15:57
看了一下基本没错的，楼主所用公式是否摘自Bolton，1980？还有，e1=lev.*q1./(0.62197+q1)+1e-10;这一条+1e ...

同有这个疑问，请问解决了吗？

logo

25行这里q1=rh1.*(0.62197.*es1./(lev-es1))/100;  
为什么不是q1=rh1.*(0.62197.*es1./(lev-0.378*es1))/100;