* Horizontal wave-activity flux derived by Takaya and Nakamura (1999, 2001)

* Used data;  Model output data of "miroc3_2_medres" for  IPCC AR4 20c3m

* Monthly-mean zonal and meridional wind (ua_A1, va_A1), geopotential height (zg_A1)

* NetCDF files have monthly-mean, multiple-level data. 

* The unit of level is [Pa], not [hPa]

* Basic state ; climatology of monthly-mean field (January)
*                   averaged from 1979 to 1998.

* Perturbation ; monthly-mean anomalies in January 1990.

* Level; 250 hPa


*-----
reinit
* u
sdfopen /h1a/IPCC-AR4/20c3m/atm/mo/ua/miroc3_2_medres/run1/ua_A1.nc
* v
sdfopen /h1a/IPCC-AR4/20c3m/atm/mo/va/miroc3_2_medres/run1/va_A1.nc
* height
sdfopen /h1a/IPCC-AR4/20c3m/atm/mo/zg/miroc3_2_medres/run1/zg_A1.nc

*  gas constant
define Ra=290
* earth radius
define a=6400000

define dlat = cdiff(lat,y)*3.14/180
define dlon = cdiff(lon,x)*3.14/180
define coslat = cos(lat*3.14/180)
define sinlat = sin(lat*3.14/180)
* Coriolis parameter
define f = 2*7.24/100000*sinlat
define f0 = 2*7.24/100000*sin(43*3.14/180)
define g=9.8

* unit [Pa]
set lev 25000

* For drawing polar projection 
set lon -5 365
set t 1

* making basic state (climatology)
define uclm = ave(ua.1,time=jan1979,time=jan1998,1yr)
define vclm = ave(va.2,time=jan1979,time=jan1998,1yr)
define zclm = ave(zg.3,time=jan1979,time=jan1998,1yr)

* anomalies
define zaa=zg.3(time=jan1990)-zclm

* QG stream function
define psiaa=g/f*zaa

* magnitude of basic state wind speed
define magU = mag(uclm,vclm)

define dpsidlon = cdiff(psiaa,x)/dlon
define ddpsidlonlon = cdiff(dpsidlon,x)/dlon

define dpsidlat = cdiff(psiaa,y)/dlat
define ddpsidlatlat = cdiff(dpsidlat,y)/dlat
define ddpsidlatlon = cdiff(dpsidlat,x)/dlon

define termxu = dpsidlon*dpsidlon-psiaa*ddpsidlonlon

*define termxv = dpsidlon*dpsidlat-ddpsidlatlon <-- !error! (2010/11/25)
define termxv = dpsidlon*dpsidlat-ddpsidlatlon*psiaa

define termyv = dpsidlat*dpsidlat-psiaa*ddpsidlatlat

* "p" is normalized by 1000hPa
define coeff1=coslat*(lev/100000)/(2*magU)
*x-component
*define px = coeff1/(a*a*coslat)*( uclm*termxu+ vclm/coslat*termxv) <-- !error!  (2010/11/25)
define px = coeff1/(a*a*coslat)*( uclm*termxu/coslat + vclm*termxv)

*y-component
define py = coeff1/(a*a)*( uclm/coslat*termxv + vclm*termyv)

set lon 0 360

set gxout contour
set cint 30
set black -0.1 0.1

* stream-function-like geopotential height
d maskout( zaa*abs(f0/f),  abs(lat)-10)

* horizontal wave-activity flux
set arrscl 0.5 20
d skip(px,6,4);maskout( py , abs(lat)-10)