哪位大神帮我看看这个假相当位温85hPa这么算对不对？谢

jiajian0512

1. 
   
2.   
   

复制代码

% 850hPa假相当位温
clear all;
clc;
grib_struct=read_grib('g:\031200',-1); % 读取NCEP再分析资料,其中-1代表全部读取
    read_grib('g:\031200','inv'); %生成变量目录表文档
   t_850hPa=grib_struct(154).fltarray;
   rh_850hPa=grib_struct(155).fltarray;
for i=1:181;
    for j=1:360;
        lon(i,j)=181-j;
        lat(i,j)=91-i;
        t(i,j)=t_850hPa((i-1)*360+j);
        rh(i,j)=rh_850hPa((i-1)*360+j);
    end
end
m_proj('equidistant','lon',[70 110],'lat',[30 70]);
prs=850;
for m=1:181;
  for n=1:360;
      es(m,n)=(6.1078*exp(17.2693882*(t(m,n)-273.16)/(t(m,n)-35.86)));
      qq(m,n)=rh(m,n)*(0.62197*es(m,n)/(prs-0.378*es(m,n)))/100;
       e(m,n)=prs*qq(m,n)/(0.62197+qq(m,n))+1e-10;
       tlcl(m,n)=55.0+2840.0/(3.5*log(t(m,n))-log(e(m,n))-4.805);
       theta(m,n)=t(m,n)*(1000/prs)^(0.2854*(1.0-0.28*qq(m,n)));
       eqt(m,n)=theta(m,n)*exp(((3376/tlcl(m,n))-2.54)*qq(m,n)*(1.0+0.81*qq(m,n)));
  end
end
  [cs,h]=m_contourf(lon,lat,eqt,'color','k');
  hold on;
  m_coast('color','k');
  m_grid;