爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 7077|回复: 17

[分享资料] 有人分析过等熵位涡么?

[复制链接]

新浪微博达人勋

发表于 2012-8-13 20:49:46 | 显示全部楼层 |阅读模式

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

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

x
有人分析过等熵位涡么?能分享下计算程序么?本人利用Grads的isen.gs老是运行出错
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2012-11-24 12:55:50 | 显示全部楼层
本帖最后由 cold_wq 于 2012-11-24 12:57 编辑


Grads自带的哈。我再分享一个最近下载的一个比较新的一个isen.gs。
  1. function isen(args)
  2. *----------------------------------------------------------------------
  3. * Bob Hart (rhart@fsu.edu) /  FSU Meteorology
  4. * Last Updated: 2/12/2012
  5. *
  6. * 2/12/12 - Updated so that it can be called as a regular function
  7. *           like cbarn.   Note that the field is now returned as
  8. *           a variable defined by the fifth argument passed. See
  9. *           updated call examples below.
  10. *
  11. * 2/26/99 - Fixed a bug that caused the script to crash on
  12. *           certain machines.  
  13. *
  14. * GrADS function to interpolate within a 3-D grid to a specified
  15. * isentropic level.  Can also be used on non-pressure level data, such
  16. * as sigma or eta-coordinate output where pressure is a function
  17. * of time and grid level.  Can be used to create isentropic PV surfaces
  18. * (examples are given at end of documentation just prior to
  19. * function.)
  20. *
  21. * Advantages:  Easy to use, no UDFs.  Disadvantages:  Can take 5-20 secs.
  22. *
  23. * Arguments:
  24. *    field = name of 3-D grid to interpolate
  25. *
  26. *    tgrid = name of 3-D grid holding temperature values (deg K) at each
  27. *            gridpoint.
  28. *
  29. *    pgrid = name of 3-D grid holding pressure values (mb) at each gridpoint
  30. *            If you are using regular pressure-level data, this should be
  31. *            set to the builtin GrADS variable 'lev'.
  32. *
  33. *    tlev  = theta-level (deg K) at which to interpolate
  34. *
  35. *    rgrid = variable name of the grid into which the result is stored
  36. *
  37. * NOTE:  Areas having tlev below bottom level or above upper level
  38. *        will be undefined in output field. Extrapolation is NOT
  39. *        performed!!
  40. *
  41. *------------------------------------------------------------------------
  42. *
  43. * EXAMPLE FUNCTION CALLS:
  44. *
  45. * Sample variables: u = u-wind in m/s
  46. *                   v = v-wind in m/s
  47. *                   w = vertical velocity
  48. *                   t = temperature in K
  49. *                   PP = pressure data in mb
  50. *
  51. * 1) Display vertical velocity field on 320K surface:
  52. *
  53. *    run isen.gs w t PP 320 omega320
  54. *    "d omega320"
  55. *
  56. * 2) Create & Display colorized streamlines on 320K surface:
  57. *
  58. *    run isen.gs u t PP 320 u320
  59. *    run isen.gs v t PP 320 v320
  60. *    "set z 1"
  61. *    "set gxout stream"
  62. *    "d u320;v320;mag(u320,v320)"
  63. *
  64. * 3) Create & display a 320K isentropic PV surface:
  65. *
  66. *    "set lev 1050 150"
  67. *    "define coriol=2*7.29e-5*sin(lat*3.1415/180)"
  68. *    "define dudy=cdiff(u,y)/(111177*cdiff(lat,y))"
  69. *    "define dvdx=cdiff(v,x)/(111177*cdiff(lon,x)*cos(lat*3.1415/180))"
  70. *    "define dt=t(z-1)*pow(1000/PP(z-1),0.286)-t(z+1)*pow(1000/PP(z+1),0.286)"
  71. *    "define dp=100*(PP(z-1)-PP(z+1))"
  72. *    "define dtdp=dt/dp"
  73. *    run isen.gs dvdx t PP 320 part1
  74. *    run isen.gs dudy t PP 320 part2
  75. *    run isen.gs dtdp t PP 320 part3
  76. *    "define pv320=-9.8*(coriol+part1-part2)*part3"
  77. *    "set z 1"
  78. *    "d pv320"
  79. *
  80. * PROBLEMS:  Send email to Bob Hart (rhart@fsu.edu)
  81. *
  82. *-----------------------------------------------------------------------
  83. *-------------------- BEGINNING OF FUNCTION ----------------------------
  84. *-----------------------------------------------------------------------

  85. * Parse input arguments

  86. field=subwrd(args,1)
  87. tgrid=subwrd(args,2)
  88. pgrid=subwrd(args,3)
  89. tlev=subwrd(args,4)
  90. rgrid=subwrd(args,5)

  91. if (rgrid = "")
  92.    say "*** NOTE:   As of March 2012, the method of calling this script has changed. ***"
  93.    say " "
  94.    say "Insufficient number of arguments given."
  95.    say " "
  96.    say "Usage:   run isen.gs field tgrid pgrid tlev rgrid"
  97.    say "    field = name of 3-D grid to interpolate"
  98.    say "    tgrid = name of 3-D grid holding temperature values (deg K) at each gridpoint."
  99.    say "    pgrid = name of 3-D grid holding pressure values (mb) at each gridpoint"
  100.    say "            If you are using regular pressure-level data, this should be set to the builtin GrADS variable 'lev'."
  101.    say "    tlev  = theta-level (deg K) at which to interpolate"
  102.    say "    rgrid = variable name of the grid into which the result is stored"
  103.    say " "
  104.    say "For example, GFS 320K zonal wind stored into the variable zonal320 would be called as:"
  105.    say "     run isen.gs ugrdprs tmpprs lev 320 zonal320"
  106.    say " "
  107.    say "Exiting."
  108.    exit
  109. endif

  110. * Get initial dimensions of dataset so that exit dimensions will be
  111. * same

  112. "q dims"
  113. rec=sublin(result,4)
  114. ztype=subwrd(rec,3)
  115. if (ztype = "fixed")
  116.    zmin=subwrd(rec,9)
  117.    zmax=zmin
  118. else
  119.    zmin=subwrd(rec,11)
  120.    zmax=subwrd(rec,13)
  121. endif

  122. * Get full z-dimensions of dataset.

  123. "q file"
  124. rec=sublin(result,5)
  125. zsize=subwrd(rec,9)

  126. * Determine spatially varying bounding pressure levels for isen surface
  127. * tabove = theta-value at level above ; tbelow = theta value at level
  128. * below for each gridpoint

  129. "set z 1 "zsize
  130. "define theta="tgrid"*pow(1000/"pgrid",0.286)"
  131. "set z 2 "zsize
  132. "define thetam="tgrid"(z-1)*pow(1000/"pgrid"(z-1),0.286)"
  133. "set z 1 "zsize-1
  134. "define thetap="tgrid"(z+1)*pow(1000/"pgrid"(z+1),0.286)"

  135. "define tabove=0.5*maskout(theta,theta-"tlev")+0.5*maskout(theta,"tlev"-thetam)"
  136. "define tbelow=0.5*maskout(theta,thetap-"tlev")+0.5*maskout(theta,"tlev"-theta)"

  137. * Isolate field values at bounding pressure levels
  138. * fabove = requested field value above isen surface
  139. * fbelow = requested field value below isen surface

  140. "define fabove=tabove*0+"field
  141. "define fbelow=tbelow*0+"field

  142. "set z 1"

  143. * Turn this 3-D grid of values (mostly undefined) into a 2-D isen layer

  144. * If more than one layer is valid (rare), take the mean of all the
  145. * valid levels. Not the best way to deal with the multi-layer issue,
  146. * but works well, rarely if ever impacts output, and is quick.
  147. * Ideally, only the upper most level would be used.  However, this
  148. * is not easily done using current GrADS intrinsic functions.

  149. "define fabove=mean(fabove,z=1,z="zsize")"
  150. "define fbelow=mean(fbelow,z=1,z="zsize")"
  151. "define tabove=mean(tabove,z=1,z="zsize")"
  152. "define tbelow=mean(tbelow,z=1,z="zsize")"

  153. * Finally, interpolate linearly in theta and create isen surface.
  154. * Linear interpolation in theta works b/c it scales as height,
  155. * or log-P, from Poisson equation for pot temp.

  156. "set z "zmin " " zmax

  157. "define slope=(fabove-fbelow)/(tabove-tbelow)"
  158. "define b=fbelow-slope*tbelow"
  159. "define "rgrid"=slope*"tlev"+b"

  160. * variable interp now holds isentropic field and its named it returned
  161. * for use by the user.

  162. say "Done.  Newly defined variable "rgrid" has "tlev"K "field"-field."

  163. return(0)
复制代码
密码修改失败请联系微信:mofangbao
回复 支持 1 反对 0

使用道具 举报

新浪微博达人勋

发表于 2012-8-14 00:54:25 | 显示全部楼层
我也用isen, 你是什么错误啊?
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2012-8-14 17:47:21 | 显示全部楼层
abd 发表于 2012-8-14 00:54
我也用isen, 你是什么错误啊?

问题已经解决了哈
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2012-8-14 20:03:25 | 显示全部楼层
楼主,能分享吗?
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2012-11-19 18:34:15 | 显示全部楼层
cold_wq 发表于 2012-8-14 17:47
问题已经解决了哈

求分享怎么用gs 处理等熵位涡  万分感谢
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2012-11-24 12:53:52 | 显示全部楼层
Heqiwei 发表于 2012-11-19 18:34
求分享怎么用gs 处理等熵位涡  万分感谢

Grads安装文件下的lib里面有个isen.gs,只是实际运用的时候,dp这一行最好这样写'define dp=100*(lev(z-1)-lev(z+1))'
,具体的例子和使用gs里面都有介绍
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2013-1-13 10:50:32 | 显示全部楼层
好东西!正好想画个等熵位涡图呢!
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2013-3-12 19:45:19 | 显示全部楼层
cold_wq 发表于 2012-11-24 12:53
Grads安装文件下的lib里面有个isen.gs,只是实际运用的时候,dp这一行最好这样写'define dp=100*(lev(z-1 ...

哈哈。。。这句很实用,我也遇到一样的问题。
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2013-3-29 15:08:17 | 显示全部楼层
下了试了能出图
密码修改失败请联系微信:mofangbao
您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

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

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

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