假相当位温脚本

紫雨微澜

分享一个假相当位温的gs的脚本，其实论坛里这类的脚本已经很多了。由于自己最近在写文章，对于公式来源等比较模糊，所以特查阅了相关参考书和文章，在此分享一下同时，还发现有些脚本还是有不少错误，当然也不完全确定这个脚本是一点错误也没有，发现错误的童鞋请指教哈~来源注释在gs附件里面，这边贴过来是乱码也附上了部分公式的参考文献。。。

1. 'reinit'
   
2. 'open **\air.ctl'
   
3. 'open **\rhum.ctl'
   
4. 'set gxout fwrite'
   
5. 'set fwrite **\thse.grd'
   
6.     tt=1
   
7.     while (tt<=365)
   
8. 'set t 'tt''
   
9.   zz=1
   
10.   while (zz<=8)
    
11. 'set z 'zz''
    
12. 'set x 1 144'
    
13. 'set  y 1 73'
    
14. 
    
15. 'define t0=air.1'
    
16. 'define rh=rhum.2'
    
17. 'define prs=lev'
    
18. 'define es=(6.1078*exp(17.2693882*(t0-273.16)/(t0-35.86)))'
    
19. *见《大气物理学》P21公式（2.2.9a）
    
20. 'define qq=rh*(0.62197*es/(prs-0.378*es))/100.'
    
21. *见《大气物理学》P20公式（2.2.5）
    
22. 'define e=prs*qq/(0.62197+qq)+1e-10'
    
23. *见《大气物理学》P20，e和r的公式自己推导一下就出来了
    
24. 'define tlcl=55.0+2840.0/(3.5*log(t0)-log(e)-4.805)'
    
25. 'define theta=t0*pow((1000/prs),(0.2854*(1.0-0.28*qq)))'
    
26. 'define eqt=theta*exp(((3376./tlcl)-2.54)*qq*(1.0+0.81*qq))'
    
27. *见《东亚副热带季风雨带建立特征及其降水性质分析》
    
28. *'define thse=eqt-273.15'
    
29. 'd eqt'
    
30. 
    
31. zz=zz+1
    
32. endwhile
    
33. tt=tt+1
    
34. endwhile  
    
35. 'disable fwrite'
    
36. 'reinit'
    
37. 
    

'reinit'<br /> 'open **\air.ctl'<br /> 'open **\rhum.ctl'<br /> 'set gxout fwrite'<br /> 'set fwrite **\thse.grd'<br />     tt=1<br />     while (tt<=365)<br /> 'set t 'tt''<br />   zz=1<br />   while (zz<=8)<br /> 'set z 'zz''<br /> 'set x 1 144'<br /> 'set  y 1 73'<br /> <br /> 'define t0=air.1'<br /> 'define rh=rhum.2'<br /> 'define prs=lev'<br /> 'define es=(6.1078*exp(17.2693882*(t0-273.16)/(t0-35.86)))'<br /> *见《大气物理学》P21公式（2.2.9a）<br /> 'define qq=rh*(0.62197*es/(prs-0.378*es))/100.'<br /> *见《大气物理学》P20公式（2.2.5）<br /> 'define e=prs*qq/(0.62197+qq)+1e-10' <br /> *见《大气物理学》P20，e和r的公式自己推导一下就出来了<br /> 'define tlcl=55.0+2840.0/(3.5*log(t0)-log(e)-4.805)' <br /> 'define theta=t0*pow((1000/prs),(0.2854*(1.0-0.28*qq)))' <br /> 'define eqt=theta*exp(((3376./tlcl)-2.54)*qq*(1.0+0.81*qq))' <br /> *见《东亚副热带季风雨带建立特征及其降水性质分析》<br /> *'define thse=eqt-273.15'<br /> 'd eqt'<br /> <br /> zz=zz+1<br /> endwhile <br /> tt=tt+1<br /> endwhile  <br /> 'disable fwrite' <br /> 'reinit'<br /> <br />

四叶草

四叶草 发表于 2015-8-18 16:26
楼主，谢谢你的分享！但其中有个公式本人感觉有点疑惑。如下：
'definee=prs*qq/(0.62197+qq)+1e-10' （1 ...

疑惑解决了！由于通常大气中的混合比和比湿都小于0.04，因此可以认为q≈r，故楼主的计算方法可取。下面是本人利用'define e=prs*qq/(0.62197+qq)+1e-10'（1）和 'define e=q*p/(0.622+0.378*q)'   （2） 出的图，图像显示几乎一样。
（1）计算的结果
（2）计算的结果

雪お轩漠

根据北大版《大气物理学》P20
我觉得程序22行的
'define e=prs*qq/(0.62197+qq)+1e-10'
应当改为
'define e=prs*p/(0.62197+qq)'
我觉得1e-10没有必要加上，我也不太明白为什么楼主加上呢，不知道我说的对不对。

pcc5201314

多谢楼主的  分享

liuzh3316

e为什么不直接es*相对湿度呢

mofangbao

遇到发帖不换行楼主可以新建立一个txt 然后把原来的文字复制进去再贴过来  已经替你改好了

紫雨微澜

mofangbao 发表于 2014-4-23 13:15
遇到发帖不换行楼主可以新建立一个txt 然后把原来的文字复制进去再贴过来  已经替你改好了

其实刚发现了，正在修改的，没想到你速度这么快。另外，怎么才能贴成这种代码形式呀~

mofangbao

紫雨微澜 发表于 2014-4-23 13:18
其实刚发现了，正在修改的，没想到你速度这么快。另外，怎么才能贴成这种代码形式呀~

发帖的时候点击代码按钮，然后选择代码为GS，在贴入文字就行了

代码高亮

加油吧，骚年！

{:eb502:}

lovechang1314

顶，楼主，正在找相关东西

wmlydia

{:eb303:}{:eb303:}终于出来了，温度的单位桑不起啊，楼主的单位是K，我的是摄氏度，公式要改一下

mai

正好要用到，谢谢！

zxc

顶，楼主，正在找相关东西，收藏了！多谢

sunground

这个是个好东西！