- 积分
- 4780
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2016-6-4
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
'reinit'
* change the directory name
* u (daily-mean field)
'sdfopen c:/1/uwnd.2017.nc'
* height (daily-mean field)
'sdfopen c:/1/hgt.2017.nc'
* gas constant
'define Ra=290'
* earth radius
'define a=6400000'
'define dlat = cdiff(lat,y)*3.1415/180'
'define dlon = cdiff(lon,x)*3.1415/180'
'define coslat = cos(lat*3.1415/180)'
'define sinlat = sin(lat*3.1415/180)'
* Coriolis parameter
'define f = 2*7.24/100000*sinlat'
'define g=9.8'
* unit [hPa]
'set lev 850'
* For drawing polar projection map
'set lon -5 365'
* Used for judgement whether the latitude is westerly or easterly
'define uclm = ave(ave(uwnd.1,lon=0,lon=360,-b),time=01sep2017,time=30sep2017)'
'define zclm = ave(ave(hgt.2,lon=0,lon=360,-b),time=01sep2017,time=30sep2017)'
* deviation from zonal-mean field
'define za=ave(hgt.2,time=01sep2017,time=30sep2017)-zclm'
* QG stream function
'define psia=g/f*za'
'define dpsidlon = cdiff(psia,x)/dlon'
'define ddpsidlonlon = cdiff(dpsidlon,x)/dlon'
'define dpsidlat = cdiff(psia,y)/dlat'
'define ddpsidlatlat = cdiff(dpsidlat,y)/dlat'
'define ddpsidlatlon = cdiff(dpsidlat,x)/dlon'
'define termx = dpsidlon*dpsidlon-psia*ddpsidlonlon'
* error, found in 2013/11/25. (Thanks to Shengping HE)
*define termy = dpsidlon*dpsidlat-psia*ddpsidlonlon
'define termy = dpsidlon*dpsidlat-psia*ddpsidlatlon'
* "p" is normalized by 1000hPa
'define coeff=(lev/1000)/(2*a*a)'
*x-component
'define px = coeff/(coslat)*(termx)'
*y-component
'define py = coeff*termy'
'set gxout contour'
*set cint 30
'set black -0.1 0.1'
* QG stream-function
'set lon 0 180'
'set lat 10 90'
'd maskout( psia, abs(lat)-10)'
* horizontal wave-activity flux
'set arrscl 0.5 10'
* maskout regions where
* (i) QG approximation is supposed to be invalid (lower latitudes)
* (ii) westerly wind speed is weak or negative
'd skip(px,2,0.5);maskout(maskout( py , abs(lat)-10),uclm-5)'
'draw title 850hpa WAF in 01-30 SEP 2017 '
'printim c:/1/DW850bo.png white'
用17年逐日的纬向风场和高度场资料,画波动通量,500hpa和200hpa的都可以出来好多纬向风矢量,但是850和925就算间隔取0.1都只能出来两行风矢量,这是为什么呢?
|
-
|