为什么我参照isen.gs的例子画位涡PV，得的每个高度上的PV都是一样的？剖面上也都直线

何小肥 · 发表于 2012-4-22 15:09:51

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

您需要登录才可以下载或查看，没有帐号？立即注册

x

现在尝试着画个等熵位涡的图，参考别的帖子里的建议，用了isen.gs，刚开始不太会就按着里面的例子来做。奇怪的是画出来的每一层的结果是一样的。。囧了

。。。为嘛？！求助！！

'reinit'
'open F:\Res\data\input\fnl\fnl_20110817_12_00.ctl'
'enable print F:\Res\draw\PV\pv315_20110817_12_00.gmf'
'set grads off'
'set grid off'
'set xlopts 1 4 0.16'
'set ylopts 1 4 0.16'
'set xlint 5'
'set ylint 5'
'set xlab on'
'set ylab on'
'set lat 20 45'
'set lon 100 130'
'set lev 1000 100'
'define u=UGRDprs'
'define v=VGRDprs'
'define PP=lev'
'define t=tmpprs'
'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,315)
'define part2='isen(dudy,t,PP,315)
'define part3='isen(dtdp,t,PP,315)
'define pv315=-9.8*(coriol+part1-part2)*part3'
'set lat 30'
'set lon 100 130'
*'set z 1'
'set cthick 6'
'cnbasemap pv315*1e6 0'
'print'
'disable print'
;
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)

'reinit' 
'open F:\Res\data\input\fnl\fnl_20110817_12_00.ctl' 
'enable print F:\Res\draw\PV\pv315_20110817_12_00.gmf' 
'set grads off' 
'set grid off' 
'set xlopts 1 4 0.16' 
'set ylopts 1 4 0.16' 
'set xlint 5' 
'set ylint 5' 
'set xlab on' 
'set ylab on' 
'set lat 20 45' 
'set lon 100 130' 
'set lev 1000 100' 
 
'define u=UGRDprs' 
'define v=VGRDprs' 
'define PP=lev' 
'define t=tmpprs' 
'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,315) 
'define part2='isen(dudy,t,PP,315) 
'define part3='isen(dtdp,t,PP,315) 
'define pv315=-9.8*(coriol+part1-part2)*part3' 
'set lat 30' 
'set lon 100 130' 
*'set z 1' 
'set cthick 6' 
'cnbasemap pv315*1e6 0' 
'print' 
'disable print' 
; 
function isen(field,tgrid,pgrid,tlev) 
*---------------------------------------------------------------------- 
* Bob Hart (<a href=\"mailto:hart@ems.psu.edu\" target=\"_blank\">hart@ems.psu.edu</a>) /  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 (<a href=\"mailto:hart@ems.psu.edu\" target=\"_blank\">hart@ems.psu.edu</a>) 
* 
*----------------------------------------------------------------------- 
*-------------------- 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)

mofangbao · 发表于 2012-4-22 21:27:25

太长了，木有用过...

何小肥 · 发表于 2012-4-23 18:53:23

mofangbao 发表于 2012-4-22 21:27

登录/注册后可看大图

太长了，木有用过...

其实37行以后的是贴的isen.gs。。。

littlelock · 发表于 2012-5-14 09:42:10

何小肥发表于 2012-4-23 18:53

登录/注册后可看大图

其实37行以后的是贴的isen.gs。。。

楼主的问题解决了没有，能用这个gs算吗？我计算出来都是恒定值。。。

pqman3 · 发表于 2012-5-14 10:56:35

你计算的是315K等熵面上的值，难道在不同高度上还有区别吗？

何小肥 · 发表于 2012-5-14 17:13:30

pqman3 发表于 2012-5-14 10:56

登录/注册后可看大图

你计算的是315K等熵面上的值，难道在不同高度上还有区别吗？

我想计算不同等压面上的值。。那我要怎么算？所有等熵面的算一遍，再在等压面上画？？用这个GS能实现么？？？

何小肥 · 发表于 2012-5-14 17:17:21

littlelock 发表于 2012-5-14 09:42

登录/注册后可看大图

楼主的问题解决了没有，能用这个gs算吗？我计算出来都是恒定值。。。

没有，最近做大尺度的，这个中小尺度的先放着。。但是听说有人用fortran做的。。只是我没有程序

littlelock · 发表于 2012-5-15 12:48:44

何小肥发表于 2012-5-14 17:17

登录/注册后可看大图

没有，最近做大尺度的，这个中小尺度的先放着。。但是听说有人用fortran做的。。只是我没有程序

我后来用那个isen.gs做出来了，只是不知道对不对

何小肥 · 发表于 2012-5-15 15:44:15

littlelock 发表于 2012-5-15 12:48

登录/注册后可看大图

我后来用那个isen.gs做出来了，只是不知道对不对

好啊。。那我找可以交流的伙伴啦~~不错不错~~有空交流下

良辰 · 发表于 2012-7-6 16:16:58

有结果了吗

		自动登录	找回密码
密码			立即注册

[分享资料] 为什么我参照isen.gs的例子画位涡PV，得的每个高度上的PV都是一样的？剖面上也都直线

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