load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
begin
;************************************************
; By YM Cheng
; Make a panel plot for the regression maps
; need regression of streamfunction, uwnd, vwnd, brtmp
;************************************************

;;=============================================
;;=========== beginning of settings ===========
;;=============================================
latRange=(/-10,35/); EOF domain
lonRange=(/-40,40/); EOF domain
PlatRange=(/-10,35/); Plot domain
PlonRange=(/-40,40/); Plot domain

diri="./regression"; directory for regression files
diro="."; directory for output figs
reg_lev=700
re_pc=1


fsf=  "sf.regres.pc"+re_pc+"."+reg_lev+"hPa.nc"
fu ="uwnd.regres.pc"+re_pc+"."+reg_lev+"hPa.nc"
fv ="vwnd.regres.pc"+re_pc+"."+reg_lev+"hPa.nc"
fbrtmp="brtmp.regres.pc"+re_pc+".nc"
;;=============================================
;;============== end of settings ==============
;;========== can leave the rest unchanged======
;;=============================================

print("================================")
print("=======plotting regression======")
print("================================")
print("fregres_sf:    "+diri+"/"+fsf)
print("fregres_u:     "+diri+"/"+fu)
print("fregres_v:     "+diri+"/"+fv)
print("fregres_brtmp: "+diri+"/"+fbrtmp)
negative=True; whether to flip the sign of regression
; negative=False
if(negative)
	sign=-1
else
	sign=1
end if

; wksType="png"
wksType="x11"


inf_sf=addfile(diri+"/"+fsf,"r")
sftmp=inf_sf->sf
sf=sftmp/100000*sign
copy_VarCoords(sftmp,sf)
delete(sftmp)


inf_u=addfile(diri+"/"+fu,"r")
utmp=inf_u->uwnd
u=utmp*sign
copy_VarMeta(utmp,u)
ndim=dimsizes(u)
delete(utmp)

inf_v=addfile(diri+"/"+fv,"r")
vtmp=inf_v->vwnd
v=vtmp*sign
copy_VarMeta(vtmp,v)
delete(vtmp)

inf_brtmp=addfile(diri+"/"+fbrtmp,"r")
brtmptmp=inf_brtmp->brtmp
brtmp=brtmptmp*sign
copy_VarMeta(brtmptmp,brtmp)
delete(brtmptmp)




mapres = True
mapres@mpMaxLatF = PlatRange(1)
mapres@mpMinLatF = PlatRange(0)
mapres@mpMaxLonF = PlonRange(1)
mapres@mpMinLonF = PlonRange(0)
mapres@gsnDraw = False
mapres@gsnFrame = False
mapres@mpOutlineOn = True
mapres@tmXBLabelFontHeightF	=	0.008
mapres@tmYLLabelFontHeightF	=	0.008
mapres@mpGeophysicalLineThicknessF = 4
mapres@mpFillOn = False
mapres@mpOutlineDrawOrder = "PostDraw"
mapres@pmTickMarkDisplayMode = "Always"
mapres@vpHeightF       = .350; set height of plot
mapres@vpWidthF        = .85	; set width of plot
mapres@vpXF            = .10; set left hand side start point of plot
mapres@vpYF		       =.75
mapres@mpGridAndLimbOn=True
mapres@mpGridLineDashPattern  = 2                 ; lat/lon lines dashed
mapres@mpGridLatSpacingF      = 5
mapres@mpGridLonSpacingF      = 10
mapres@gsnMaximize         = True


cnres=True
cnres@gsnFrame =       False
cnres@gsnDraw  =       False
cnres@cnLineLabelPlacementMode =       "Constant"
cnres@cnFillOn = False
cnres@cnLinesOn        =       True                  ; turn off contour line
cnres@tmXBMode		=	"Automatic"
cnres@cnInfoLabelOn  = True
cnres@vpXF            = .10; set left hand side start point of plot
cnres@vpYF		       =.75
cnres@gsnContourNegLineDashPattern = 1
cnres@gsnContourZeroLineThicknessF=0
cnres@cnLineThicknessF =       1
cnres@cnLevelSelectionMode ="ManualLevels"
cnres@gsnRightString	=""
cnres@gsnLeftString	=""
cnres@cnInfoLabelOn = False


res                       = True                 ; plot mods desired
res@vcRefMagnitudeF       = 1                 ; define vector ref mag
res@vcRefLengthF          = 0.035                ; define length of vec ref
res@vcGlyphStyle          = "LineArrow"        ; turn on curley vectors
res@gsnRightString=""
res@vcRefAnnoOrthogonalPosF = -.183
res@vcRefAnnoParallelPosF   = .995
res@vcRefAnnoPerimOn          = False
res@vcRefAnnoString2On        = False
res@vcLineArrowHeadMinSizeF=0.000
res@vcLineArrowHeadMaxSizeF=0.005
res@vcLineArrowThicknessF   = 1.5               ; change vector thickness
res@gsnDraw               = False                ; don't draw yet
res@gsnFrame              = False                ; don't advance frame yet
res@gsnLeftString         = " "                   ; turn off strings

shres=True
shres@gsnFrame =       False
shres@gsnDraw  =       False
shres@lbLabelBarOn         = False       ; turn off individual lb's
shres@cnLineLabelPlacementMode =       "Constant"
shres@cnFillOn = True
shres@cnLinesOn        =       False                  ; turn off contour line
shres@tmXBMode		=	"Automatic"
shres@vpXF            = .10; set left hand side start point of plot
shres@vpYF		       =.75
shres@cnLevelSelectionMode = "ExplicitLevels"

shres@gsnLeftString=""
shres@gsnRightString=""
shres@lbLabelAutoStride=False
shres@gsnAddCyclic = True
shres@tiMainFontHeightF = 0.020
shres@gsnStringFontHeightF = 0.015
shres@cnInfoLabelOn=False



wksName=diro+"/regression-v"
wks	=	gsn_open_wks(wksType,wksName)
gsn_define_colormap(wks,"BlWhRe")
shres@cnLevels=(/-5,-4.5,-4,-3.5,-3,-2.5,-2,-1.5,-1,-0.5,0.5,1,1.5,2,2.5,3,3.5,4,4.5,5/)


shres@gsnLeftString="Domain: ("\
+latRange(0)+","+latRange(1)+")/("+lonRange(0)+","+lonRange(1)+"),PC"+re_pc

shres@gsnRightString=reg_lev+" hPa sf (10~S~5~N~ m~S~2~N~s~S~-1~N~); wind vector (ms~S~-1~N~); Tb (K)"
cnres@cnMaxLevelValF=20; sf
cnres@cnMinLevelValF=-20
cnres@cnLevelSpacingF=2


i=40; Day -2
lag_day=u&lag(i)
centstring="Day "+sprintf("%5.2f",lag_day)
print("Lag "+lag_day+" now")
shres@tiMainString=centstring
plot1	:=	gsn_csm_map(wks,mapres)
plot2 := gsn_csm_contour(wks,sf(i,:,:),cnres)
plot3   :=gsn_csm_vector(wks,u(i,:,:),v(i,:,:),res)
plot4	:=gsn_csm_contour(wks,brtmp(i,:,:),shres)
overlay(plot1,plot4)
overlay(plot1,plot2)
overlay(plot1,plot3)
shres@gsnRightString=""
shres@gsnLeftString=""

i=48; Day 0
lag_day=u&lag(i)
centstring="Day "+sprintf("%5.2f",lag_day)
print("Lag "+lag_day+" now")
shres@tiMainString=centstring
plot5	:=	gsn_csm_map(wks,mapres)
plot6 := gsn_csm_contour(wks,sf(i,:,:),cnres)
plot7 :=gsn_csm_vector(wks,u(i,:,:),v(i,:,:),res)
plot8	:=gsn_csm_contour(wks,brtmp(i,:,:),shres)
overlay(plot5,plot8)
overlay(plot5,plot7)
overlay(plot5,plot6)


i=56; Day 2
lag_day=u&lag(i)
centstring="Day "+sprintf("%5.2f",lag_day)
print("Lag "+lag_day+" now")
shres@tiMainString=centstring
plot9	:=	gsn_csm_map(wks,mapres)
plot10 := gsn_csm_contour(wks,sf(i,:,:),cnres)
plot11   :=gsn_csm_vector(wks,u(i,:,:),v(i,:,:),res)
plot12	:=gsn_csm_contour(wks,brtmp(i,:,:),shres)
overlay(plot9,plot12)
overlay(plot9,plot11)
overlay(plot9,plot10)



pnl = True
pnl@amJust = "TopLeft"
pnl@gsnPanelFigureStringsFontHeightF = 0.010
pnl@gsnPanelFigureStringsPerimOn = False
pnl@amOrthogonalPosF = -0.496
pnl@amParallelPosF = -0.496
pnl@gsnPanelLabelBar           = True    ; Turn on common labelbar
pnl@lbLabelAutoStride          = True    ; Spacing of lbar labels.
pnl@lbLabelFontHeightF = 0.008
pnl@gsnMaximize=True
pnl@lbLabelAutoStride=False
pnl@lbLabelStride=1
gsn_panel(wks,(/plot1,plot5,plot9/),(/3,1/),pnl); frame(wks)

end