;This script divides EPz by lev, but not scale EPy and EPz by sqrt(pho)
;plot pdf figure for models and differences in one file.
;filtered EP flux to the power of 1/3
;*******************************************************
; Code written by Joe Barsugli joseph.barsugli@noaa.gov
; from the NOAA/ESRL PSD
; adapted by Cathy Smith cathy.smoth@noaa.gov
; last checked Nov 2009.
; modified by Joe Barsugli to add contours of EP-Flux divergence June 2010
; modified by Joe Barsugli to redo scaling of arrows in the vertical June 2010
;********************************************************
;
; +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"      	; functions required to
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"           	; plot.  include before
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"         	; begin
; +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
;
begin

;Name1="spcamctl"
;Name2="WA1Jan_FREE"
;Case1="spcamctl"
;Case2="IDEAL_WA1Jan_FREE"
Name1="CONTROL_noQBO"
Name2="dmpdz2e-3_noQBO"
Case1="WACCMSC_CONTROL_noQBO"
Case2="WACCMSC_dmpdz2e-3_noQBO"
;Case1="IDEAL_WA1Jan_FREE"
;Case2="IDEAL_WA1Jan_EQ1M40_5"
;Name1="WA1Jan_FREE"
;Name2="WA1Jan_EQ1M40_5"
;Case1="IDEAL_WA1Jan_FREE"
;Case2="IDEAL_WA1Jan_EQ1M40lon60180_5"
;Name1="WA1Jan_FREE"
;Name2="WA1Jan_EQ1M40lon60180_5"
;Case1="IDEAL_WA1Jan2D_FREE"
;Case2="IDEAL_WA1Jan2D_EQ1M40_5"
;Name1="WA1Jan2D_FREE"
;Name2="WA1Jan2D_EQ1M40_5"
Name1="CONTROL_noQBO"
Name2="dmpdz2e-3_noQBO"
Case1="WACCMSC_CONTROL_noQBO"
Case2="WACCMSC_dmpdz2e-3_noQBO"
;Case1="ERA40"
;Case2="spcam_nlev51"
;Name1="ERA40"
;Name2="SP-Strat"
;Case1="IDEAL_SP1JAN_FREE_ENS"
;Case2="IDEAL_SP1JAN_EQ1M40"
;Name1="SP1JAN_FREE_ENS"
;Name2="SP1JAN_EQ1M40"
;Case1="IDEAL_WA1Jan_FREE_BLK_1000200_1545"
;Case2="IDEAL_WA1Jan_EQ1M40_BLK_1000200_1545"
;Name1="WA1Jan_FREE_BLK_1000200_1545"
;Name2="WA1Jan_EQ1M40_BLK_1000200_1545"
Case1="IDEAL_WA1Jan35_FREE"
Case2="IDEAL_WA1Jan35_EQ1M40_10"
Name1="WA1Jan35_FREE"
Name2="WA1Jan35_EQ1M40_10"
;Case1="IDEAL_WA1Jan35_FREE"
;Case2="IDEAL_WA1Jan35_EQ1M40lon60180_5"
;Name1="WA1Jan35_FREE"
;Name2="WA1Jan35_EQ1M40lon60180_5"
;Case1="IDEAL_WA1Jan10m35_EQ1M40_2"
;Case2="IDEAL_WA1Jan10m35_EQ1M40_2"
;Name1="WA1Jan10m35_EQ1M40_2"
;Name2="WA1Jan10m35_EQ1M40_2"
Case1="IDEAL_WA1Jan10m35_FREE"
Case2="IDEAL_WA1Jan10m35_EQ1M40_2"
Name1="WA1Jan10m35_FREE"
Name2="WA1Jan10m35_EQ1M40_2"
;Case1="IDEAL_WA1Jan35_FREE"
;Case2="IDEAL_WA1Jan35_EQ1M40_10"
;Name1="WA1Jan35_FREE"
;Name2="WA1Jan35_EQ1M40_10"

Case_tmp=getenv("Case2")
if (.not. ismissing(Case_tmp)) then
	Case2=Case_tmp
end if
Name_=str_split(Case2,"_")
Name2=str_join(Name_(1:),"_")
print(Name2)

IDEAL=True
season="DJF"

sf=1
scale_by_sqrt_p=0
showtop=4
showbot=500
showbot=900

showdiffdivbot=500.
showdiffepbot=500.
;showdiffdivbot=600.
;showdiffepbot=600.
;showdiffdivbot=990.
;showdiffepbot=990.
largerscale=5.
nocontour=False
weaksource=False
scalefilter=3

lefton=True
righton=True

;ptype="waccm"
;level0=(/0,0/)
;vinthkeeplev0=level0-(/0,0/)
;nlev=66 - max(level0)

ptype="waccm"
level0=(/31,31/)
vinthkeeplev0=level0-(/31,31/)
nlev=66 - max(level0)

;ptype="cam30"
;level0=(/0,0/)
;vinthkeeplev0=level0-(/0,0/)
;nlev=30 - max(level0)


;filters=(/"","stationary","filter","filter","filter","filter","filter","filter"/)
;freq1s=(/"","","20","20","1","1","100","100"/)
;freq2s=(/"","","100","100","20","20","10000","10000"/)
;addstationarys=(/"","","sta","","sta","","sta",""/)
;scale_1_3s=(/False,False,False,False,False,False,False,False/)
;scale_1_3s=(/True,True,True,True,True,True,True,True/)


filters=(/"","stationary","filter","filter","filter","filter","filter","filter","filter","filter","filter","filter","filter","filter","filter","filter","filter","filter","filter","filter","filter","filter","other","other"/)
freq1s=(/"","","35","35","1","1","1","1","45","45","40","40","40","40","40","40","40","40","40","0","0","1","40","35"/)
freq2s=(/"","","45","45","35","35","7","7","1000","1000","40","40","40","40","40","40","40","40","40","10000","10000","15","40","45"/)
addstationarys=(/"","","sta","","sta","","sta","","sta","","","","","","","","","","","","","","minus","minus"/)
;scale_1_3s=(/False,False,False,False,False,False,False,False,False,False,False,False,False,False,False,False,False,False,False,False,False,False,False,False/)
scale_1_3s=(/True,True,True,True,True,True,True,True,True,True,True,True,True,True,True,True,True,True,True,True,True,True,True,True,True/)
ffts=(/"kf","kf","kf","kf","kf","kf","kf","kf","kf","kf","single","single","single","single","single","single","single","single","single","k1","k2","kf","single","kf"/)
appendEPnames=(/"","","","","","","","","","","","construct","construct_k2","construct_k1","stationary","stationary_k2","stationary_k1","k1","k2","","","","",""/)
appendEPnameouts=(/"","","","","","","","","","","","construct","construct_k2","construct_k1","stationary","stationary_k2","stationary_k1","k1","k2","k1","k2","","",""/)
;"","construct","construct_k2","stationary","stationary_k2" just used in 4040-single filtered files

;do iplt=0,dimsizes(filters)-2
;do iplt=dimsizes(filters)-1,dimsizes(filters)-1
;do iplt=1,1 ; stationary
do iplt=0,0
;do iplt=3,3 ; 35-45
;do iplt=8,8
;do iplt=10,10 ;40-40
filter=filters(iplt)
freq1=freq1s(iplt)
freq2=freq2s(iplt)
addstationary=addstationarys(iplt)
scale_1_3=scale_1_3s(iplt)
fft=ffts(iplt)
; "";"fft";"kf";"single","k1","k2"
appendEPname=appendEPnames(iplt)
appendEPnameout=appendEPnameouts(iplt)

paperdir="paper2/"
;paperdir=Name1+"_"+Name2+"/"
diri="./T_U_V_uv/"
diriscratch="/glade/scratch/wykang/output/"
diro="./EPflux/"
Names=(/Name1,Name2/)
append=""
if (scale_1_3) then
	append="_p13"
end if
if (nocontour) then
	append=append+"_nocontour"
end if

if (filter.eq."") then
freq1=""
freq2=""
; if no filter, set freq1 and freq2 to ""
end if

totalcase=2-1

show_accel=0

strat=""
if(showbot.le. 100) then
  strat="strat"
end if

f_1=addfile(diri+"yz_climatology_"+Name1+season+".nc","r")
f_2=addfile(diri+"yz_climatology_"+Name2+season+".nc","r")
if (fileexists(diri+"vinth_UVTuvvt_"+Name1+season+".nc")) then
	f1=addfile(diri+"vinth_UVTuvvt_"+Name1+season+".nc","r")
	f2=addfile(diri+"vinth_UVTuvvt_"+Name2+season+".nc","r")
else
	f1=addfile(diri+"vinth_climatology_"+Name1+season+".nc","r")
	f2=addfile(diri+"vinth_climatology_"+Name2+season+".nc","r")
end if


if (filter .eq. "filter" .or. filter.eq."other") then
;	feddy_1=addfile(diriscratch+Case1+"/"+Case1+"_filter_"+freq1+freq2+season+fft+".nc","r")
;	feddy_2=addfile(diriscratch+Case2+"/"+Case2+"_filter_"+freq1+freq2+season+fft+".nc","r")
       feddy_1=addfile(diri+Case1+"_filter_"+freq1+freq2+season+fft+".nc","r")
       feddy_2=addfile(diri+Case2+"_filter_"+freq1+freq2+season+fft+".nc","r")
print(diri+Case1+"_filter_"+freq1+freq2+season+fft+".nc")
		vtname="v_t_"+appendEPname
        vt_1=(/feddy_1->$vtname$(vinthkeeplev0(0):,:)/)
        vt_2=(/feddy_2->$vtname$(vinthkeeplev0(1):,:)/)
		uvname="u_v_"+appendEPname
        uv_1=(/feddy_1->$uvname$(vinthkeeplev0(0):,:)/)
        uv_2=(/feddy_2->$uvname$(vinthkeeplev0(1):,:)/)
end if
if (filter.eq."") then
	feddy_1=f_1
	feddy_2=f_2
	vt_1=feddy_1->v_t_
	vt_2=feddy_2->v_t_
	uv_1=feddy_1->u_v_
	uv_2=feddy_2->u_v_
end if
if (filter.eq."other") then
	feddy_1=f_1
	feddy_2=f_2
	print(dimsizes(feddy_1->v_t_))
	print(dimsizes(vt_1))
	vt_1=feddy_1->v_t_-vt_1

	vt_2=feddy_2->v_t_-vt_2
	uv_1=feddy_1->u_v_-uv_1
	uv_2=feddy_2->u_v_-uv_2
end if

;if (freq1.eq."20".and.freq2.eq."100") then
;
;       ;feddy_1=addfile(diri+Case1+"_filter_10010000"+season+fft+".nc","r")
;       ;feddy_2=addfile(diri+Case2+"_filter_10010000"+season+fft+".nc","r")
;
;	   ;vt_1=vt_1+feddy_1->v_t_
;	   ;vt_2=vt_2+feddy_2->v_t_
;	   ;uv_1=uv_1+feddy_1->u_v_
;	   ;uv_2=uv_2+feddy_2->u_v_
;
;
;       feddy_1=addfile(diri+Case1+"_filter_135"+season+fft+".nc","r")
;       feddy_2=addfile(diri+Case2+"_filter_135"+season+fft+".nc","r")
;		vt_1=f_1->v_t_-doubletofloat(dim_avg((f1->T)*(f1->V))-(f_1->T)*(f_1->V))-feddy_1->v_t_
;		vt_2=f_2->v_t_-doubletofloat(dim_avg((f2->T)*(f2->V))-(f_2->T)*(f_2->V))-feddy_2->v_t_
;		uv_1=f_1->u_v_-doubletofloat(dim_avg((f1->U)*(f1->V))-(f_1->U)*(f_1->V))-feddy_1->u_v_
;		uv_2=f_2->u_v_-doubletofloat(dim_avg((f2->U)*(f2->V))-(f_2->U)*(f_2->V))-feddy_2->u_v_
;       feddy_1=addfile(diri+Case1+"_filter_3545"+season+fft+".nc","r")
;       feddy_2=addfile(diri+Case2+"_filter_3545"+season+fft+".nc","r")
;	   vt_1=vt_1-feddy_1->v_t_
;	   t_2=vt_2-feddy_2->v_t_
;	   uv_1=uv_1-feddy_1->u_v_
;	   uv_2=uv_2-feddy_2->u_v_
;
;end if
;
if (filter.eq."stationary") then
	vt_1=doubletofloat(dim_avg((f1->T)*(f1->V))-(f_1->T)*(f_1->V))
        vt_2=doubletofloat(dim_avg((f2->T)*(f2->V))-(f_2->T)*(f_2->V))
	uv_1=doubletofloat(dim_avg((f1->U)*(f1->V))-(f_1->U)*(f_1->V))
        uv_2=doubletofloat(dim_avg((f2->U)*(f2->V))-(f_2->U)*(f_2->V))
end if

if (addstationary.eq."sta") then
	vt_1=vt_1+doubletofloat(dim_avg((f1->T)*(f1->V))-(f_1->T)*(f_1->V))
	vt_2=vt_2+doubletofloat(dim_avg((f2->T)*(f2->V))-(f_2->T)*(f_2->V))
	uv_1=uv_1+doubletofloat(dim_avg((f1->U)*(f1->V))-(f_1->U)*(f_1->V))
	uv_2=uv_2+doubletofloat(dim_avg((f2->U)*(f2->V))-(f_2->U)*(f_2->V))
end if
if (addstationary.eq."minus") then
	vt_1=vt_1 - doubletofloat(dim_avg((f1->T)*(f1->V))-(f_1->T)*(f_1->V))
	vt_2=vt_2 - doubletofloat(dim_avg((f2->T)*(f2->V))-(f_2->T)*(f_2->V))
	uv_1=uv_1 - doubletofloat(dim_avg((f1->U)*(f1->V))-(f_1->U)*(f_1->V))
	uv_2=uv_2 - doubletofloat(dim_avg((f2->U)*(f2->V))-(f_2->U)*(f_2->V))
end if

lat=f_1->lat

T_bar_1=f_1->T
T_bar_2=f_2->T

;T_bar_CESM1=f_CESM1->T
;T_bar_CESM4=f_CESM4->T

U_bar_1=f_1->U
U_bar_2=f_2->U

uv_str=0
st="uv_"+Name1
uv_str@$st$=uv_1
st="uv_"+Name2
uv_str@$st$=uv_2
;uv_str@uv_CESM1=uv_CESM1
;uv_str@uv_CESM4=uv_CESM4

vt_str=0
st="vt_"+Name1
vt_str@$st$=vt_1
st="vt_"+Name2
vt_str@$st$=vt_2

T_bar_str=0
st="T_bar_"+Name1
T_bar_str@$st$=T_bar_1
st="T_bar_"+Name2
T_bar_str@$st$=T_bar_2
;T_bar_str@T_bar_CESM1=T_bar_CESM1
;T_bar_str@T_bar_CESM4=T_bar_CESM4

U_bar_str=0
st="U_bar_"+Name1
U_bar_str@$st$=U_bar_1
st="U_bar_"+Name2
U_bar_str@$st$=U_bar_2

Fj=(/uv_1(0:29,0:29),uv_2(0:29,0:29)/)
Fk=(/vt_1(0:29,0:29),vt_2(0:29,0:29)/)
U=(/U_bar_1(0:29,0:29),U_bar_2(0:29,0:29)/)
DIV=U
plotEP=new(3,graphic)

;!!!!!!!!!!!start loop!!!!!!!!!!!!!!
do n = 0,totalcase,1
st="T_bar_"+Names(n)
T_bar=T_bar_str@$st$
st="U_bar_"+Names(n)
U_bar = U_bar_str@$st$

incoordinate=addfile("/glade/scratch/wykang/analysis/vertical_coordinate_"+ptype+".nc","r")
pres = incoordinate->lev(level0(n):)
pres@units="mb"
pres@long_name="Pressure(mbar)"
lev=pres
;H=10000
H=1
theta_bar = T_bar
theta_bar=doubletofloat(T_bar*(conform(T_bar,pres,0)/1000.)^(-0.286))

theta_p = theta_bar
theta_p_temp = center_finite_diff_n(theta_bar,log(pres),False,0,0)
theta_p_temp = theta_p_temp/conform(theta_p_temp,pres*100.,0)
theta_p_temp = theta_p_temp/H
theta_p=doubletofloat(theta_p_temp)

U_p=center_finite_diff_n(U_bar,log(pres),False,0,0)/H/conform(theta_p_temp,pres*100.,0)

st="vt_"+Names(n)


v_t_bar=vt_str@$st$
v_theta_bar=v_t_bar
v_theta_bar=doubletofloat(v_t_bar*(conform(v_t_bar,pres,0)/1000.)^(-0.286))
st="uv_"+Names(n)
u_v_bar=uv_str@$st$


a = 6.37122e06 ;  radius of the earth
;a=1;
PI=3.14159265358979
phi = lat*PI/180.0     ; latitude in radians
acphi=a*cos(phi)       ; 
asphi=a*sin(phi)       ; a* sin latitude for use in calculating the divergence.
omega = 7.2921e-5      ;
f = 2*omega*sin(phi)   ; coriolis parameter

latfac=acphi*cos(phi)/a   ; scale factor includes extra cos(phi) for graphical display of arrows (see Edmon et al, 1980)
;U_y=center_finite_diff_n(U_bar*conform(U_bar,cos(phi),1),conform(U_bar,a*phi,1),False,0,1)
U_y=center_finite_diff_n(U_bar,conform(U_bar,a*phi,1),False,0,1)
;F_j = (-u_v_bar)*conform(u_v_bar,latfac,1)
;F_k = (conform(v_theta_bar,f*cos(phi),1))*v_theta_bar/theta_p*conform(u_v_bar,100*pres,0)

theta_p(20:,45:51)=where(theta_p(20:,45:51).le.2e-4,2e-4,theta_p(20:,45:51))
F_j = (U_p*v_theta_bar/theta_p-u_v_bar)*conform(u_v_bar,latfac,1)
;F_k = (conform(v_theta_bar,f*cos(phi),1)+U_y)*v_theta_bar/theta_p*conform(u_v_bar,100*pres,0)
F_k = (conform(v_theta_bar,f*cos(phi),1)+U_y*conform(U_bar,cos(phi),1))*v_theta_bar/theta_p
;F_k = (conform(v_theta_bar,f*cos(phi),1))*v_theta_bar/theta_p
;here Fk is density weighted!
copy_VarMeta(u_v_bar,F_j)
copy_VarMeta(v_theta_bar,F_k)
F_j@long_name="EP flux phi"
F_k@long_name="EP flux p"
F_j!0="lev"
F_j!1="lat"
F_j&lev=pres
F_j&lat=lat
F_k!0="lev"
F_k!1="lat"
F_k&lev=pres
F_k&lat=lat
U_bar!0="lev"
U_bar!1="lat"
U_bar&lev=pres
U_bar&lat=lat
F_k_=F_k
F_k_=-F_k
;**************************************************************
;save
;**************************************************************

system("rm "+diro+"yz_EP_"+filter+"_"+freq1+freq2+"_"+season+"_"+Names(n)+appendEPnameout+".nc")
fout=addfile(diro+"yz_EP_"+filter+"_"+freq1+freq2+"_"+season+"_"+Names(n)+appendEPnameout+".nc","c")
fout->F_j=F_j
fout->F_k=F_k
delete(fout)

;**************************************************************
;divergence and interpolate to standard grid
;**************************************************************
; take derivative w.r.t latitude  using  1/[a cos(phi)]  d/dphi [cos(phi)*X] =  d/d [asin(phi)] F_j
; note that F_j already has the extra  factor of cos(phi)
  EPdiv1  = center_finite_diff(F_j,asphi,False,0)
; take derivative w.r.t pressure  (in Pascals, hence the factor of 100.0)
;  Need to re-order the coordinates to pass to center_finite_diff 
  EPdiv2temp  = center_finite_diff(F_k(lat|:,lev|:),lev*100.0,False,0);
  EPdiv2temp!0="lat" 
  EPdiv2temp!1="lev"
  EPdiv2 = (/EPdiv2temp(lev|:,lat|:) /)  ; And then re-order the coordiantes again  (NCL awkwardness! will be fixed in NCL 5.2)
  EPdiv = (EPdiv1 + EPdiv2)*86400    ;   Add these together to get the divergence of F
  copy_VarMeta(F_k,EPdiv)

; interpolate to a regular grid in log(pressure) ; This dimension re-ordering due to NCL awkwardness
; show  every other latitude -- excluding the poles. 
lev_int = 10^fspan(3,0,30)   ; interpolation targets
lev_int@units="hpa"
lat_int = fspan(-90,90,30)
linlog=-2                       ; Option to int2p that gives log-interpolation with no extrapolation
F_k_int3 = int2p(lev,F_k(lat|:,lev|:),lev_int,linlog)
F_j_int3 = int2p(lev,F_j(lat|:,lev|:),lev_int,linlog)
EPdiv_int3 = int2p(lev,EPdiv(lat|:,lev|:),lev_int,linlog)
U_bar_int3 = int2p(lev,U_bar(lat|:,lev|:),lev_int,linlog)
F_k_int3!0="lat"
F_k_int3!1="lev"
F_j_int3!0="lat"
F_j_int3!1="lev"
EPdiv_int3!0="lat"
EPdiv_int3!1="lev"
U_bar_int3!0="lat"
U_bar_int3!1="lev"
F_k_int2 = int2p(lat,F_k_int3(lev|:,lat|:),lat_int,1)
F_j_int2 = int2p(lat,F_j_int3(lev|:,lat|:),lat_int,1)
EPdiv_int2 = int2p(lat,EPdiv_int3(lev|:,lat|:),lat_int,1)
U_bar_int2 = int2p(lat,U_bar_int3(lev|:,lat|:),lat_int,1)
F_k_int2!0="lev_int"
F_k_int2&lev_int=lev_int
F_k_int2!1="lat_int"
F_k_int2&lat_int=lat_int
copy_VarMeta(F_k_int2,F_j_int2)
copy_VarMeta(F_k_int2,EPdiv_int2)
copy_VarMeta(F_k_int2,U_bar_int2)

; re-order
F_k_int = F_k_int2(lev_int|:,lat_int|:)
F_j_int = F_j_int2(lev_int|:,lat_int|:)
EPdiv_int = EPdiv_int2(lev_int|:,lat_int|:)
U_bar_int = U_bar_int2(lev_int|:,lat_int|:)

if (show_accel .eq. 1) then
  dudt_int = 86400.0*EPdiv_int/conform(EPdiv_int,a*cos(phi(1:71:2)),1) ; 
  copy_VarMeta(EPdiv_int,dudt_int)
else
  dudt_int = 0.0*EPdiv_int
end if

; Scale the vectors for display 
; First scale according to Edmon et al. for pressure coordinates
;      (even though I am using log-p display -- not entirely consistent as the arrows may "look" divergent when they are not, but better visibility in practice)
F_k_int = F_k_int   
F_j_int = F_j_int/a*30

; Next scale by the relative ranges of the two axes of the plot(  3.14 radians by 10^5 Pa)
F_k_int = F_k_int/conform(F_k_int,lev_int*100,0)/(10^0.1-1)
F_j_int = F_j_int/PI

; Optionaly scale by sqrt(presure)
if (scale_by_sqrt_p .eq. 1)
print("scale!")
rhofac = sqrt(1000/lev_int)
F_k_int = F_k_int*conform(F_k_int,rhofac,0)
F_j_int = F_j_int*conform(F_j_int,rhofac,0)
end if


; Scale by a magnification factor above 100 mb. 
strat1 = new(30,float)
strat1 = (/ 0., 1., 1., 1., 1., 1., 1.,1.,1., 1., 1., sf,sf,sf,sf,sf,sf,sf,sf,sf,sf,sf,sf,sf,sf,sf, sf, sf, sf, sf/)
stratmask=conform(F_k_int,strat1,0)
F_k_int = F_k_int*stratmask
F_j_int = F_j_int*stratmask

;clean up not-shown area
showtoplev=ind(lev_int .le. showtop)
if (.not.ismissing(showtoplev(0))) then
F_j_int(showtoplev,:)=0
F_k_int(showtoplev,:)=0
EPdiv_int(showtoplev,:)=0
end if
showbotlev=ind(lev_int .ge. showbot)
if (.not.ismissing(showbotlev(0))) then
F_j_int(showbotlev,:)=0
F_k_int(showbotlev,:)=0
EPdiv_int(showbotlev,:)=0
end if

; Hide the 1000 mb lev for the divergence and acceleration
dudt_int@_FillValue = -999.0
dudt_int(0,:)=dudt_int@_FillValue
EPdiv_int@_FillValue = -999.0
EPdiv_int(0,:)=EPdiv_int@_FillValue

;remove too large data
F_j_int=where(abs(F_j_int).ge.0.0003.or.abs(F_k_int).ge.0.0003,0,F_j_int )
F_k_int=where(abs(F_j_int).ge.0.0003.or.abs(F_k_int).ge.0.0003,0,F_k_int )
EPdiv_int=where(abs(EPdiv_int).ge.10., 0,EPdiv_int)

Fj(n,:,:)=doubletofloat(F_j_int)
Fk(n,:,:)=doubletofloat(F_k_int)

;scale by taking 1/3 power
if (scale_1_3) then
F_k_int = doubletofloat(((F_k_int^2+F_j_int^2)*10^10)^(-1./3.)*F_k_int)
F_j_int = doubletofloat(((F_k_int^2+F_j_int^2)*10^10)^(-1./3.)*F_j_int)
F_j_int = where(isnan_ieee(F_j_int),0,F_j_int)
F_k_int = where(isnan_ieee(F_k_int),0,F_k_int)
end if

;************************************************
; Create Plot 
;************************************************
vectitle =   Names(n)+" "+season+"   "+freq1+"-"+freq2

; create vector plot resources for pressure-lev grid  (not used for plotting in this version)
  res_vec_int = True
;  res_vec_int@tiMainString          = vectitle                     ; plot title
  res_vec_int@tiMainFontHeightF     = 0.025
  res_vec_int@tiXAxisFontHeightF    = 0.025
  res_vec_int@tiYAxisFontHeightF    = 0.015
  res_vec_int@tiXAxisString         = "latitude"           ; x-axis label
if (lefton) then
  res_vec_int@tiYAxisString         = "pressure (mb)"         ; y-axis label
  res_vec_int@tmYLMode        = "Explicit"
  res_vec_int@tmYLValues      = (/showtop,10,30,50,100,300,500,800/)
  res_vec_int@tmYLLabels	=(/showtop,10,30,50,100,300,500,800/)
else
  res_vec_int@tiYAxisString         = ""
  res_vec_int@tmYLOn               = False 
  res_vec_int@tmYLLabelsOn         = False  
end if

  res_vec_int@trYReverse            = True                 ; reverse y-axis
  res_vec_int@gsnYAxisIrregular2Log = True                 ; set y-axis to log scale
  res_vec_int@trYMinF    =  showtop
  res_vec_int@trYMaxF    =  showbot
  res_vec_int@trXMinF = -90.
  res_vec_int@vfYArray              = lev_int
  res_vec_int@vfXArray              = lat_int


  res_vec_int@vcGlyphStyle		= "LineArrow"
  res_vec_int@vcMonoLineArrowColor	= True
  res_vec_int@vcRefMagnitudeF       = 0.0001
  res_vec_int@vcRefLengthF          = 0.04               

if (.not.scale_1_3) then
	res_vec_int@vcRefMagnitudeF       =res_vec_int@vcRefMagnitudeF*1.
end if

;if (filter .eq. "filter") then
;  res_vec_int@vcRefMagnitudeF=res_vec_int@vcRefMagnitudeF/scalefilter
;end if

  res_vec_int@vcFillArrowFillColor  =False
  res_vec_int@vcFillArrowEdgeColor = False

  res_vec_int@pmLabelBarDisplayMode = "NoCreate"          ; Turn on a label bar.
;  res_vec_int@pmLabelBarWidthF      = 0.08              ; make it thinner
;  res_vec_int@lbPerimOn             = False             ; no box around it
  res_vec_int@vpWidthF              = 0.60
  res_vec_int@vpHeightF             = 0.35
;  res_vec_int@lbLabelBarOn          = False                     ; turn off label bar
;  res_vec_int@gsnMaximize           = True

 res_vec_int@gsnDraw                 = False
 res_vec_int@gsnFrame                = False  
 res_vec_int@vcRefAnnoOn             = False
; res_vec_int@lbOrientation              ="Vertical"

; Create contour plot resources for interpolated grid
  res_con_int                       = True
  res_con_int@trYMinF    	    =  showtop
  res_con_int@trYMaxF    	    =  showbot
  res_con_int@sfYArray              = lev_int                ; use pressure for y axis
  res_con_int@sfXArray              = lat_int             ; use lat for x axis
  res_con_int@cnFillOn              = True     ; turn on color fill
  res_con_int@trYReverse            = True                 ; reverse y-axis
  res_con_int@gsnYAxisIrregular2Log = True                 ; set y-axis to log scale

 res_con_int@gsnDraw                = False
 res_con_int@gsnFrame               = False

 res_con_int@cnLinesOn       	    = False
 res_con_int@gsnContourZeroLineThicknessF = 2.0
 res_con_int@gsnContourPosLineDashPattern = 0
 res_con_int@gsnContourNegLineDashPattern = 2
 res_con_int@cnSmoothingOn          = True
 res_con_int@cnLineLabelsOn         = False
 res_con_int@gsnContourLineThicknessesScale  = 1

 res_con_int@cnLevelSelectionMode = "ManualLevels"
 res_con_int@cnMinLevelValF	   =-10.
 res_con_int@cnMaxLevelValF		=10.
 res_con_int@cnLevelSpacingF		=1.0

if (filter.eq."filter") then
 res_con_int@cnMinLevelValF=res_con_int@cnMinLevelValF/scalefilter
 res_con_int@cnMaxLevelValF=res_con_int@cnMaxLevelValF/scalefilter
 res_con_int@cnLevelSpacingF=res_con_int@cnLevelSpacingF/scalefilter
end if

if (.not.righton) then
 res_con_int@pmLabelBarDisplayMode = "Never"
else
 res_con_int@pmLabelBarDisplayMode = "Always"
end if

 ;res_con_int@lbLabelBarOn          = righton
 res_con_int@lbPerimOn              	= False
; res_con_int@lbOrientation		="Horizontal"
; res_con_int@lbLabelPosition		="Bottom"
 res_con_int@cnInfoLabelOn = False
res_con_int@pmLabelBarWidthF        = 0.1

if (lefton) then
  res_con_int@tiYAxisString         = "pressure (mb)"         ; y-axis label
  res_con_int@tmYLMode        = "Explicit"
  res_con_int@tmYLValues      = (/showtop,10,30,50,100,300,500,800/)
  res_con_int@tmYLLabels        =(/showtop,10,30,50,100,300,500,800/)
else
  res_con_int@tiYAxisString         = ""
  res_con_int@tmYLOn               = False
  res_con_int@tmYLLabelsOn         = False
end if

if (n.eq.0) then
  wks  = gsn_open_wks("pdf",diro+paperdir+Names(0)+"_"+Names(1)+"_yz_EP_"+filter+freq1+freq2+addstationary+"_"+season+strat+appendEPname+append)
  gsn_define_colormap(wks,"BlWhRe") 
end if
  plotvec = gsn_vector(wks,F_j_int,F_k_int,res_vec_int)     ; creates plot
  ;printMinMax(F_j_int,True)
  ;printMinMax(F_k_int,True)
  plotEP(n) = gsn_contour(wks,EPdiv_int,res_con_int)     ; creates plot for div(F)
  overlay(plotEP(n),plotvec)
 draw(plotEP(n))
 frame(wks)

;prepare for next loop
;Fj(n,:,:)=doubletofloat(F_j_int)
;Fk(n,:,:)=doubletofloat(F_k_int)
U(n,:,:)=doubletofloat(U_bar_int)
DIV(n,:,:)=doubletofloat(EPdiv_int)
delete(T_bar)
delete(U_bar)
delete(pres)
delete(lev)
delete(theta_bar)
delete(theta_p)
delete(theta_p_temp)
delete(U_p)
delete(v_t_bar)
delete(u_v_bar)
delete(v_theta_bar)
delete(U_y)
delete(F_j)
delete(F_k)
delete(F_k_)
delete(EPdiv1)
delete(EPdiv2)
delete(EPdiv2temp)
delete(EPdiv)
end do
;===========================================================
; Compare the two cases
;===========================================================
dFj=F_j_int
dFk=F_k_int
dU=U_bar_int
dDIV=EPdiv_int

dFj=(/rm_single_dims(Fj(1,:,:)-Fj(0,:,:))/)
dFk=(/rm_single_dims(Fk(1,:,:)-Fk(0,:,:))/)
dU=(/rm_single_dims(U(1,:,:)-U(0,:,:))/)
dDIV=(/rm_single_dims(DIV(1,:,:)-DIV(0,:,:))/)

;scale by taking 1/3 power
if (scale_1_3) then
dFk = doubletofloat(((dFk^2+dFj^2)*10^10)^(-1./3.)*dFk)
dFj = doubletofloat(((dFk^2+dFj^2)*10^10)^(-1./3.)*dFj)
dFj = where(isnan_ieee(dFj),0,dFj)
dFk = where(isnan_ieee(dFk),0,dFk)
end if

copy_VarMeta(F_j_int,dFj)
copy_VarMeta(F_k_int,dFk)
copy_VarMeta(U_bar_int,dU)
copy_VarMeta(EPdiv_int,dDIV)

F_j_int=(/dFj/)
F_k_int=(/dFk/)
EPdiv_int=(/dDIV/)

showdiffdivbotlev=ind(lev_int.ge. showdiffdivbot)
dDIV(showdiffdivbotlev,:)=0
showdiffepbotlev=ind(lev_int.ge. showdiffepbot)
dFj(showdiffepbotlev,:)=0
dFk(showdiffepbotlev,:)=0
;remove too large data
;dFj=where(abs(dFj).ge.0.0001.or.abs(dFk).ge.0.0001,0,dFj )
;dFk=where(abs(dFj).ge.0.0001.or.abs(dFk).ge.0.0001,0,dFk )
;dDIV=where(abs(dDIV).ge.3., 0,dDIV)


;=============================plot difference===============================

print("plot difference")
  res_vec_int@vcRefMagnitudeF       = 0.00002
if (.not.scale_1_3) then
	res_vec_int@vcRefMagnitudeF       =res_vec_int@vcRefMagnitudeF*1.
end if
;if (filter .eq. "filter") then
;  res_vec_int@vcRefMagnitudeF=res_vec_int@vcRefMagnitudeF/scalefilter
;end if
if (IDEAL.and.weaksource) then
tropolev=ind(lev_int.gt.60.)
dFj(tropolev,:)=dFj(tropolev,:)/20.
dFk(tropolev,:)=dFk(tropolev,:)/20.
end if

  plotvec = gsn_vector(wks,dFj,dFk,res_vec_int)     ; creates plot

  res_con_int@cnLevelSelectionMode = "ManualLevels"

if (IDEAL.and.weaksource) then
tropolev=ind(lev_int.gt.60.)
dDIV(tropolev,:)=dDIV(tropolev,:)/10.
end if

  res_con_int@cnLevelSpacingF=0.02*largerscale
  res_con_int@cnMaxLevelValF=  0.19*largerscale
  res_con_int@cnMinLevelValF=  -0.19*largerscale

;if (filter .eq. "filter") then
;  res_con_int@cnLevelSpacingF=res_con_int@cnLevelSpacingF/scalefilter
;  res_con_int@cnMaxLevelValF=res_con_int@cnMaxLevelValF/scalefilter
;  res_con_int@cnMinLevelValF=res_con_int@cnMinLevelValF/scalefilter
;  res_con_int@cnLevelSpacingF=0.01
;  res_con_int@cnMaxLevelValF=  0.095
;  res_con_int@cnMinLevelValF=  -0.095
;end if

if (nocontour) then
  plotEP(2) = gsn_contour(wks,dDIV*0.,res_con_int)     ; creates plot for div(F)
else
  plotEP(2) = gsn_contour(wks,dDIV,res_con_int)     ; creates plot for div(F)
end if
  overlay(plotEP(2),plotvec)
 draw(plotEP(2))
 frame(wks)

  panel_res=True

  panel_res@gsnMaximize         = True         ; large format
  panel_res@gsnPanelBottom      = 0.05         ; add some space at bottom
  panel_res@amJust   = "TopLeft"
  panel_res@gsnPanelFigureStringsFontHeightF = 0.004
  panel_res@gsnPanelFigureStrings = (/"~F22~"+Name1,"~F22~"+Name2,"~F33~D~F22~EP"/)
  panel_res@gsnPanelLabelBar    = False         ; add common colorbar
;  panel_res@lbLabelBarOn	= True
;  panel_res@lbLabelAutoStride   = True         ; auto stride on labels
;  panel_res@lbLabelFontHeightF  = 0.01
  gsn_panel(wks,plotEP,(/2,2/),panel_res)


end do;iplt
end