波动通量的计算

董伟是我

'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都只能出来两行风矢量，这是为什么呢？

river

你不是用了maskout函数吗，其他的可能都是不符合条件的，所以就没画出来吧

董伟是我

river 发表于 2018-4-17 18:43
你不是用了maskout函数吗，其他的可能都是不符合条件的，所以就没画出来吧

对对，非常感谢，就是这个问题