等熵位涡 求指导错误

Joucly

'reinit'
'open e:\home\fnl_20130629_00_00.ctl'
'set lon 90 140'
'set lat 20 60'
'set lev 1000 100'
'set t 1'
'define tt=TMPprs'
'define u=ugrdprs'
'define v=vgrdprs'

'define coriol=2*7.29e-5*sin(lat*3.1415/180)'
'define dudy=cdiff(u,y)/(111177*cdiff(lat,y))'
'define dvdx=cdiff(v,x)/(111177*cdiff(lon,x)*cos(lat*3.1415/180))'
'define dt=tt(z-1)*pow(1000/lev(z-1),0.286)-tt(z+1)*pow(1000/lev(z+1),0.286)'
'define dp=100*(lev(z-1)-lev(z+1))'
'define dtdp=dt/dp'
'define part1='isen(dvdx,tt,lev,320)
'define part2='isen(dudy,tt,lev,320)
'define part3='isen(dtdp,tt,lev,320)
'define pv320=-9.8*(coriol+part1-part2)*part3'
'set z 1'
'set lon 100 130'
'set lat 30 50'
'd pv320*1000000'

'printim e:\630\pv\dsang\1.jpg x1024 y768 white'

*后面都是自带的gs
function isen(field,tgrid,pgrid,tlev)
*----------------------------------------------------------------------
* Bob Hart (hart@ems.psu.edu) /  PSU Meteorology
* 2/26/1999
*
* 2/26/99 - Fixed a bug that caused the script to crash on
*           certain machines.  
*
* GrADS function to interpolate within a 3-D grid to a specified
* isentropic level.  Can also be used on non-pressure level data, such
* as sigma or eta-coordinate output where pressure is a function
* of time and grid level.  Can be used to create isentropic PV surfaces
* (examples are given at end of documentation just prior to
* function.)
*
* Advantages:  Easy to use, no UDFs.  Disadvantages:  Can take 5-20 secs.
*
* Arguments:
*    field = name of 3-D grid to interpolate
*
*    tgrid = name of 3-D grid holding temperature values (deg K) at each
*            gridpoint.
*
*    pgrid = name of 3-D grid holding pressure values (mb) at each gridpoint
*            If you are using regular pressure-level data, this should be
*            set to the builtin GrADS variable 'lev'.
*
*    tlev  = theta-level (deg K) at which to interpolate
*
* Function returns:  defined grid interp holding interpolated values
*
* NOTE:  YOU NEED TO INCLUDE A COPY OF THIS FUNCTION IN YOUR SCRIPT
*
* NOTE:  Areas having tlev below bottom level or above upper level
*        will be undefined in output field. Extrapolation is NOT
*        performed!!
*
*------------------------------------------------------------------------
*
* EXAMPLE FUNCTION CALLS:
*
* Sample variables: u = u-wind in m/s
*                   v = v-wind in m/s
*                   w = vertical velocity
*                   t = temperature in K
*                  PP = pressure data in mb
*
* 1) Display vertical velocity field on 320K surface:
*
*    "d "isen(w,t,PP,320)
*
* 2) Create & Display colorized streamlines on 320K surface:
*
*    "define u320="isen(u,t,PP,320)
*    "define v320="isen(v,t,PP,320)
*    "set z 1"
*    "set gxout stream"
*    "d u320;v320;mag(u320,v320)"
*
* 3) Create & display a 320K isentropic PV surface:
*
*    "set lev 1050 150"
*    "define coriol=2*7.29e-5*sin(lat*3.1415/180)"
*    "define dudy=cdiff(u,y)/(111177*cdiff(lat,y))"
*    "define dvdx=cdiff(v,x)/(111177*cdiff(lon,x)*cos(lat*3.1415/180))"
*    "define dt=t(z-1)*pow(1000/PP(z-1),0.286)-t(z+1)*pow(1000/PP(z+1),0.286)"
*    "define dp=100*(PP(z-1)-PP(z+1))"
*    "define dtdp=dt/dp"
*    "define part1="isen(dvdx,t,PP,320)
*    "define part2="isen(dudy,t,PP,320)
*    "define part3="isen(dtdp,t,PP,320)
*    "define pv320=-9.8*(coriol+part1-part2)*part3"
*    "set z 1"
*    "d pv320"
*
* PROBLEMS:  Send email to Bob Hart (hart@ems.psu.edu)
*
*-----------------------------------------------------------------------
*-------------------- BEGINNING OF FUNCTION ----------------------------
*-----------------------------------------------------------------------

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

"q dims"
rec=sublin(result,4)
ztype=subwrd(rec,3)
if (ztype = "fixed")
   zmin=subwrd(rec,9)
   zmax=zmin
else
   zmin=subwrd(rec,11)
   zmax=subwrd(rec,13)
endif

* Get full z-dimensions of dataset.

"q file"
rec=sublin(result,5)
zsize=subwrd(rec,9)

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

"set z 1 "zsize
"define theta="tgrid"*pow(1000/"pgrid",0.286)"
"set z 2 "zsize
"define thetam="tgrid"(z-1)*pow(1000/"pgrid"(z-1),0.286)"
"set z 1 "zsize-1
"define thetap="tgrid"(z+1)*pow(1000/"pgrid"(z+1),0.286)"

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

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

"define fabove=tabove*0+"field
"define fbelow=tbelow*0+"field

"set z 1"

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

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

"define fabove=mean(fabove,z=1,z="zsize")"
"define fbelow=mean(fbelow,z=1,z="zsize")"
"define tabove=mean(tabove,z=1,z="zsize")"
"define tbelow=mean(tbelow,z=1,z="zsize")"

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

"set z "zmin " " zmax

"define slope=(fabove-fbelow)/(tabove-tbelow)"
"define b=fbelow-slope*tbelow"
"define interp=slope*"tlev"+b"

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

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

return(interp)




运行出来的错误

平流层的萝卜

我也在算位涡，不过还没用isen插值到等熵面上。楼主不妨看看未插值前的850hPa的位涡如何，看看位涡算的是否正确

平流层的萝卜

还有，你是先算各个分量，把它们插值之后再一起算的位涡。有没有试过直接把位涡算出来，再插值到等熵面上呢？

平流层的萝卜

咱们用的资料正好是一个，这是我算的广义湿位涡——http://bbs.06climate.com/forum.php?mod=viewthread&tid=17828，但是不是用isen.gs里边的算法算的。一起交流下~

Joucly

平流层的萝卜 发表于 2013-10-22 21:06

                               
登录/注册后可看大图

还有，你是先算各个分量，把它们插值之后再一起算的位涡。有没有试过直接把位涡算出来，再插值到等熵面上呢 ...

谢谢啊····

刚刚自己搞定了···

我正要用广义位涡就看见你的帖子······ 现在正在用你的gs试呢

Joucly

平流层的萝卜 发表于 2013-10-22 21:06

                               
登录/注册后可看大图

还有，你是先算各个分量，把它们插值之后再一起算的位涡。有没有试过直接把位涡算出来，再插值到等熵面上呢 ...

楼主问你个问题····  你帖子里提到的吴国雄的文章  能告诉我是哪篇么？ 我看了很多广义位涡的文章 都没有你贴出来的这么详尽  大多数直接给出的公式推导结果

平流层的萝卜

Joucly 发表于 2013-10-22 21:23

                               
登录/注册后可看大图

楼主问你个问题····  你帖子里提到的吴国雄的文章  能告诉我是哪篇么？ 我看了很多广义位涡的文章 都 ...

涡旋发展和移动的动力和热力问题Ⅰ-_i_PV-Q__i_观点.pdf。从这个里边截图的，详细的推导我没推过，水平有限。

平流层的萝卜

Joucly 发表于 2013-10-22 21:18

                               
登录/注册后可看大图

谢谢啊····

刚刚自己搞定了···

哦？是什么问题呢？

Joucly

平流层的萝卜 发表于 2013-10-22 21:27

                               
登录/注册后可看大图

哦？是什么问题呢？

set lev 的问题···  grads里面不能经度  维度  高度  都是变化的

Joucly

Joucly 发表于 2013-10-22 21:23

                               
登录/注册后可看大图

楼主问你个问题····  你帖子里提到的吴国雄的文章  能告诉我是哪篇么？ 我看了很多广义位涡的文章 都 ...

谢谢啊···  我看看呢··