load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"

begin

; for use of high resolution map database, execute
; 'setenv NCARG_RANGS /network/rit/lab/atm419lab/zhou/rangs' (on snow)
; 'setenv NCARG_RANGS /atm419lab/zhou/rangs' (on ash, reed)
; before running the script

c="./cmap_uv200.txt"
cmap=asciiread(c,(/83,3/),"float")
cmap=cmap/255

data=addfile("./UVZ_composite_sub.nc","r")

;U=data->U
;U!0="case"
;U!1="lat"
;U!2="lon"

;V=data->V
;V!0="case"
;V!1="lat"
;V!2="lon"

Z=data->Z
Z!0="case"
Z!1="lat"
Z!2="lon"

;MAG=sqrt(U^2+V^2)

Z&case=ispan(0,52,1)
Z&lat=fspan(-20,20,57)
Z&lon=fspan(-20,20,57)

;copy_VarMeta(U,V)
;copy_VarMeta(U,Z)
;copy_VarMeta(U,MAG)

dsizes=dimsizes(Z)
;print(dsizes)
dcase=dsizes(0)
dlat=dsizes(1)
dlon=dsizes(2)

nEOF=3
optEOF=True
optEOF@jopt=0
;U_EOF=eofunc_n_Wrap(U,nEOF,optEOF,0)
;pcvar_U=U_EOF@pcvar
;V_EOF=eofunc_n_Wrap(V,nEOF,optEOF,0)
;pcvar_V=V_EOF@pcvar
Z_EOF=eofunc_n_Wrap(Z,nEOF,optEOF,0)
pcvar_Z=Z_EOF@pcvar

;printVarSummary(U_EOF)
;print(U_EOF@pcvar)

;printVarSummary(V_EOF)
;print(V_EOF@pcvar)

printVarSummary(Z_EOF)
print(Z_EOF@pcvar)

;MAG_EOF=sqrt(U_EOF^2+V_EOF^2)
;copy_VarMeta(U_EOF,MAG_EOF)

;printVarSummary(MAG_EOF)
;print(MAG_EOF@pcvar)

;f=5*10^-5
;g=9.81

U_EOF=new((/nEOF,dlat,dlon/),"float")
V_EOF=new((/nEOF,dlat,dlon/),"float")

do n=0,nEOF-1,1
  do j=1,dlat-2,1
    do i=1,dlon-2,1

      U_EOF(n,j,i)=-10*(Z_EOF(n,j+1,i)-Z_EOF(n,j-1,i))/2.0
      V_EOF(n,j,i)=10*(Z_EOF(n,j,i+1)-Z_EOF(n,j,i-1))/2.0

    end do
  end do
end do

MAG_EOF=sqrt(U_EOF^2+V_EOF^2)

copy_VarCoords(Z_EOF,U_EOF)
copy_VarCoords(Z_EOF,V_EOF)
copy_VarCoords(Z_EOF,MAG_EOF)

do n=0,nEOF-1,1

print("EOF "+(n+1)+" U range: "+min(U_EOF(n,:,:))+" "+max(U_EOF(n,:,:)))
print("EOF "+(n+1)+" V range: "+min(V_EOF(n,:,:))+" "+max(V_EOF(n,:,:)))
print("EOF "+(n+1)+" Z range: "+min(Z_EOF(n,:,:))+" "+max(Z_EOF(n,:,:)))
print("EOF "+(n+1)+" MAG range: "+min(MAG_EOF(n,:,:))+" "+max(MAG_EOF(n,:,:)))

wks=gsn_open_wks("png","eof_zgrad_cov_"+(n+1))
gsn_define_colormap(wks,cmap)

res=True

res@gsnAddCyclic=False
res@gsnFrame=False
res@gsnDraw=False
res@gsnMaximize=False
res@gsnStringFont=21
res@gsnLeftString="200hPa Geopotential Height"
res@gsnLeftStringOrthogonalPosF=0.01
res@gsnLeftStringFontHeightF=0.015
res@gsnRightString="EOF "+(n+1)+" Z "+sprintf("%3.1f",pcvar_Z(n))+"%" ;+\
                                ;" V "+sprintf("%3.1f",pcvar_V(n))+"%"+\
				;" Z "+sprintf("%3.1f",pcvar_Z(n))+"%"
;res@gsnRightString="composite of 92 cases"
;res@gsnRightString="+"+sprinti("%0.2i",counter)+"h"
;res@gsnRightString="+"+sprintf("%4.1f",counter*0.5)+"h"
res@gsnRightStringOrthogonalPosF=0.01
res@gsnRightStringFontHeightF=0.015
res@gsnScalarContour=True

res@vpXF=0.15
res@vpYF=0.95
res@vpHeightF=0.7
res@vpWidthF=0.7

res@trYMinF=-20
res@trYMaxF=20
res@trXMinF=-20
res@trXMaxF=20

res@vcRefMagnitudeF=0.06
res@vcRefLengthF=0.035
res@vcGlyphStyle="CurlyVector"
res@vcLineArrowThicknessF=3.0
res@vcRefAnnoOrthogonalPosF=-0.089
res@vcRefAnnoParallelPosF=1.0
res@vcMinDistanceF=0.03
res@vcMonoLineArrowColor=True
res@vcLineArrowColor=1
res@vcRefAnnoArrowUseVecColor=False
res@vcRefAnnoArrowLineColor=1
res@vcLineArrowHeadMinSizeF=0.003
res@vcLineArrowHeadMaxSizeF=0.008

res@vcRefAnnoString2On=False

;res@stLineColor=0
;res@stArrowLengthF=0.004
;res@stMinArrowSpacingF=0.004
;res@stMinDistanceF=0.1
;res@stLineThicknessF=2.0
;res@stStepSizeF=0.01

res@cnFillOn=True
res@cnLevelSelectionMode="ExplicitLevels"
res@cnLevels:=fspan(0,0.07,71)
res@cnFillOpacityF=1
res@cnLinesOn=False

res@lbBoxSeparatorLinesOn=False
res@lbLabelStride=10

res@pmLabelBarHeightF=0.1

;res@mpDataResolution="FinestResolution"
;res@mpDataBaseVersion="HighRes"
;res@mpLimitMode="LatLon"
;;res@mpCenterLatF=centlat
;;res@mpCenterLonF=centlon
;res@mpMinLatF=centlat-20
;res@mpMaxLatF=centlat+20
;res@mpMinLonF=centlon-20
;res@mpMaxLonF=centlon+20
;res@mpMonoFillColor=False
;res@mpFillColors=(/-1,-1,-1/)
;res@mpOutlineOn=True

res@tmXTOn                = False
res@tmYROn                = False

res@tmXBMode              = "Explicit"
res@tmXBValues            = ispan(-15,15,5)
res@tmXBLabels            = res@tmXBValues
;res@tmXBLabels            = (/"100~S~o~N~E","110~S~o~N~","120~S~o~N~E",\
;                              "130~S~o~N~","140~S~o~N~E","150~S~o~N~",\
;                              "160~S~o~N~E","170~S~o~N~","180~S~o~N~E"/)
res@tmXBMajorLengthF      = 0.01
res@tmXBMajorOutwardLengthF=0.0
res@tmXBMajorThicknessF   = 3.0
res@tmXBMinorOn           = False
res@tmXBMinorPerMajor     = 4
res@tmXBMinorLengthF      = 0.004
res@tmXBMinorOutwardLengthF=0.0
res@tmXBMinorThicknessF   = 2.0
res@tmXBLabelFont         = 21
res@tmXBLabelFontHeightF  = 0.015
res@tmXBLabelDeltaF       = -0.6

res@tmYLMode              = "Explicit"
res@tmYLValues            = ispan(-15,15,5)
res@tmYLLabels            = res@tmYLValues
;res@tmYLLabels            = (/"0~S~o~N~","10~S~o~N~N","20~S~o~N~",\
;                              "30~S~o~N~N","40~S~o~N~","50~S~o~N~N"/)
res@tmYLMajorLengthF      = 0.01
res@tmYLMajorOutwardLengthF=0.0
res@tmYLMajorThicknessF   = 3.0
res@tmYLMinorOn           = False
res@tmYLMinorPerMajor     = 4
res@tmYLMinorLengthF      = 0.004
res@tmYLMinorOutwardLengthF=0.0
res@tmYLMinorThicknessF   = 2.0
res@tmYLLabelFont         = 21
res@tmYLLabelFontHeightF  = 0.015
res@tmYLLabelDeltaF       = -0.3

plot1=gsn_csm_contour(wks,MAG_EOF(n,:,:),res)

res@cnFillOn=False
res@cnLinesOn=True
res@cnInfoLabelOn=False

if (n.eq.0) then
  Z_EOF(n,:,:)=-Z_EOF(n,:,:)
  U_EOF(n,:,:)=-U_EOF(n,:,:)
  V_EOF(n,:,:)=-V_EOF(n,:,:)
end if

if (n.eq.1) then
  Z_EOF(n,:,:)=-Z_EOF(n,:,:)
  U_EOF(n,:,:)=-U_EOF(n,:,:)
  V_EOF(n,:,:)=-V_EOF(n,:,:)
end if

if (n.eq.2) then
  Z_EOF(n,:,:)=-Z_EOF(n,:,:)
  U_EOF(n,:,:)=-U_EOF(n,:,:)
  V_EOF(n,:,:)=-V_EOF(n,:,:)
end if

;if (n.eq.2) then
;  Z_EOF(n,:,:)=-Z_EOF(n,:,:)
;end if


zmin=tointeger(min(Z_EOF(n,:,:))*1000)
zmax=tointeger(max(Z_EOF(n,:,:))*1000)



res@cnLevels:=ispan(floattoint(zmin),floattoint(zmax),5)/1000.0
res@cnLineLabelPlacementMode="constant"
res@cnLineLabelInterval=2
res@cnLineThicknessF=3.0

plot2=gsn_csm_contour(wks,smth9_Wrap(Z_EOF(n,:,:),0.5,0.25,False),res)

;if (n.eq.0) then
; ; U_EOF(n,:,:)=-U_EOF(n,:,:)
;  V_EOF(n,:,:)=-V_EOF(n,:,:)
;end if


;if (n.eq.1) then
;  U_EOF(n,:,:)=-U_EOF(n,:,:)
;  V_EOF(n,:,:)=-V_EOF(n,:,:)
;end if

;if (n.eq.2) then
;  ;U_EOF(n,:,:)=-U_EOF(n,:,:)
;  V_EOF(n,:,:)=-V_EOF(n,:,:)
;end if


plot3=gsn_csm_vector(wks,U_EOF(n,:,:),V_EOF(n,:,:),res)

draw(plot1)
draw(plot2)
draw(plot3)

frame(wks)

end do

end