爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 46374|回复: 35

[脚本编辑] 关于露点温度计算的探讨(已知干球温度和相对湿度)

[复制链接]

新浪微博达人勋

发表于 2017-3-25 23:08:36 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
本帖最后由 river 于 2021-9-11 00:15 编辑

   (2021年9月11日更新新版图片,由于之前出图以后没有处理,导致上传论坛以后版面很难看。这次利用了一页多图,这样图片排列整齐有利于对比,顺便把所用的脚本进行更新,希望可以帮到大家,也希引起大家进一步思考。)
    家园里有很多关于露点温度计算的帖子,但是使用的公式都没有明确说明,而且有同一个人发的帖子,但是计算公式的系数上有差别,不清楚是什么问题。比如这个帖子http://bbs.06climate.com/forum.php?mod=viewthread&tid=38232,提出了问题,但是最终也没有得出结论;另外还有这个帖子http://bbs.06climate.com/forum.p ... 0%B6%D4%CA%AA%B6%C8,其中三楼 兰北 同学给出了一个露点温度计算公式,但是不清楚出处,而且计算结果有问题,先不考虑;六楼 西风江南 给出了VC环境下的程序,应该是利用了迭代法,但是我能力有限,暂时还没想出怎么在GrADS脚本里实现。

   我自己又仔细看了看南信大高庆九老师的《天气学诊断分析》,里面有详细的关于温湿参量的计算公式,于是我就借助这个课件上的公式自己进行了相关公式的推导,最终也得到了一个露点计算的方法。(其实课件上也是说要通过迭代法计算露点温度,但是我试了半天算出来最大的能到35度,而且很大一个范围,所以应该还是有问题,后面再试试);我还在一个Python画斜温图的脚本里找到了一个计算露点的公式,放在GrADS脚本里一起和 @笨笨不笨 @传说中的谁 一起作比较。所以就是目前一共有四种公式来计算露点温度。

现在我把以上几种计算方法计算得到的同一天同一时次的500/700/850/1000hPa的露点温度结果放在下面,供大家讨论:

500.png
700.png
850.png
1000.png




从图中不难看出,其实四种方法计算出来的露点大同小异,整体趋势是一致的,个别细节有差异。而且越往高空差别越大,是因为越往高空温度越低,而低于冰点的露点计算本来就是目前的一大难点,所以不用太纠结。至于使用哪一种公式,如果是基本高于冰点的情况下的露点计算使用哪一种都一样,低于冰点的情况就需要稍微斟酌一下。个人感觉还是 传说中的谁 的公式计算出来的露点温度分布更均匀一些,不像前两个计算出来的个别地方特别密集!

最后就是期待有大神把迭代法计算露点的方法写进GrADS的脚本里!

下面是我使用的脚本,有问题请大家及时指出!
  1. 'reinit'
  2. 'open g:/ncep/2015/201504.ctl'

  3. fig ='Td1 Td2 Td3 Td4'

  4. lev.1=500;lev.2=700;lev.3=850;lev.4=1000;

  5. tt=1
  6. while(tt<=4)
  7. 'set t 3'

  8. 'set lev 'lev.tt''

  9. *-------已知温度、相对湿度、气压计算比湿--------
  10. 'define TC=TMPprs-273.16' ;*某高度层的摄氏温度℃
  11. 'define TK=TMPprs' ;*某高度层的开氏温度K
  12. 'define RH=RHprs' ;*某高度层的相对湿度Relative humidity%
  13. 'define PRS=lev' ;*获得某层高度的气压hPa
  14. *-------水面-------
  15. if(TC>=-15)
  16. 'AA=17.2693882'
  17. 'BB=35.86'
  18. endif
  19. *-------冰面-------       
  20. if(TC<=-40)
  21. 'AA=21.8745584'
  22. 'BB=7.66'
  23. endif          
  24. *-----冰水混合面----
  25. if(-40<=TC<=-15)        
  26. 'define AA=21.8745584+(17.2693882-21.8745584)/(-15-(-40.))*(TC-(-40))'
  27. 'define BB=7.66+(35.86-7.66)/(-15-(-40.))*(TC-(-40))'
  28. endif
  29. *-----求饱和水汽压Tetens经验公式hPa--------
  30. 'define ES=6.1078*exp(AA*TC/(TK-BB))'
  31. ******水汽压hPa******
  32. 'define e=ES*RH/100+1e-10'

  33. *if(e<ES)
  34. *'define TC=TC-0.05'
  35. *'define ES=6.1078*exp(AA*TC/(TK-BB))'
  36. *else
  37. *'define td=Tc'
  38. *endif

  39. 'define td1=( (273.16-BB)*log(e)+1.80957*BB-494.3012 )/(AA-log(e)+1.80957)'

  40. ******Use Bolton's (1980, MWR, p1047) formulae to find tdew****
  41. 'define ratio=log(e/6.112)'
  42. 'define Td2=(((17.67-ratio)*273.16+243.5*ratio)/(17.67-ratio))-273.16'
  43.    
  44. 'define td3=TC-((14.55+0.114*TC)*(1-0.01*RH) + pow((2.5+0.007*TC)*(1-0.01*RH),3) + (15.9+0.117*TC)*pow((1-0.01*RH),14))'

  45. 'define td4=TC-((14.55+0.114*TC)*(1-0.01*RH) + pow((2.5+0.007*TC)*(1-0.01*RH),3) + (15.9+0.37*TC)*pow((1-0.01*RH),14))'

  46. *******图形输出******
  47. i=1
  48. while(i<=4)

  49. 'set lon 85 125'
  50. 'set lat 30 55'
  51. 'set mpdset cnworld'
  52. 'set map 15 1 6'
  53. *'set mproj scaled'
  54. *'set parea 1 10 1 8'
  55. 'set grads off'
  56. 'set grid off'

  57. 'set ylopts 1 5 0.15'
  58. 'set xlopts 1 5 0.15'
  59. 'set annot 1 5'

  60. 'set vpage 0 11 0 8.5'
  61. 'set grads off'

  62. 'subplot_y 2 2 'i''

  63. 'set ylint 5'
  64. 'set xlint 5'

  65. *'set csmooth on'

  66. 'set gxout shaded'
  67. *'set cthick 7'
  68. *'set cstyle 3'
  69. *'set ccolor 2'
  70. 'set cint 5'
  71. 'run G:\GrADS-NCL-colormakerV1.3\output\23colorsBlueToRed.gs'
  72. 'd smth9(Td'i')'
  73. *'cbarn 0.65 1 7.2 5'

  74. 'set gxout contour'
  75. 'set cthick 8'
  76. 'set cint 5'
  77. 'set ccolor 1'
  78. *'set cstyle 1'
  79. *'set clopts 1'
  80. *'set clskip 1'
  81. *'set clab masked'
  82. *'set clab forced'
  83. 'd smth9(Td'i')'

  84. * draw the map
  85. 'set mpdraw on'
  86. 'set mpdset cnhimap'
  87. 'set map 15 1 6'
  88. 'draw map'
  89. *'huanghe.gs'
  90. *'rivers.gs'

  91. * 'draw_figstr2 TL 1.0 Fig.'i''
  92. figi=subwrd(fig,i)
  93. 'draw_figstr2 TL 1.0 ('figi')'

  94. i=i+1
  95. endwhile
  96. 'set string 1 c 7 0'
  97. 'set strsiz 0.2 0.2'
  98. 'draw string 5.5 8 'lev.tt'hPa Dewpoint'
  99. *'cbar_interp 1 1 1 9.7 4.3'
  100. 'gxprint G:\ncep\2015\pic\'lev.tt'.png x2400 y1854 white'
  101. 'c'
  102. tt=tt+1
  103. endwhile
  104. ;

复制代码









来自群组: 龙王山皇家气象学院

评分

参与人数 3金钱 +45 贡献 +20 收起 理由
我想要台时光机 + 5 很给力!
mofangbao + 20 + 20
四叶草 + 20 赞一个!

查看全部评分

密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2017-3-25 23:11:43 | 显示全部楼层
本帖最后由 river 于 2017-3-26 20:34 编辑

不知道图片怎么能弄成一排四个那样的呢@兰溪之水@言深深 @传说中的谁 @mofangbao
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2017-3-25 23:47:05 | 显示全部楼层
谢谢楼主!先收着,还在学习阶段
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2017-3-27 09:02:15 | 显示全部楼层
谢谢楼主分享!平时计算都是直接用文献里的公式,没深究
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2017-3-27 11:39:27 | 显示全部楼层
一页多图啊,清风发过帖子的
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2017-3-27 11:56:39 | 显示全部楼层
传说中的谁 发表于 2017-3-27 11:39
一页多图啊,清风发过帖子的

我是说论坛发帖子这个功能只能按照原图片大小,不能进行其他排版操作好像······
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2017-3-28 19:13:10 | 显示全部楼层
先收藏再学习!
密码修改失败请联系微信:mofangbao

新浪微博达人勋

0
早起挑战累计收入
发表于 2017-3-29 12:51:04 | 显示全部楼层
river 发表于 2017-3-27 11:56
我是说论坛发帖子这个功能只能按照原图片大小,不能进行其他排版操作好像······

你的原图太大了。。。
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2017-3-31 09:29:10 | 显示全部楼层
mofangbao 发表于 2017-3-29 12:51
你的原图太大了。。。

嗯,确实是。看后面有空的话弄小一点,或者一页多图可能更好
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2017-4-1 08:49:32 | 显示全部楼层
mark
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表