;
;    avganomplot_moistflx.ncl
;
;
;    This program takes precipitation data in from generated netCDF files which
;    were created by timegrab.ncl, and processed by catprecip.csh and pstats.ncl.
;    The data from these files are then plotted for the appropriate 
;    region of North/Central America and the Pacific Ocean plus Gulf Stream. 
;
;    This program plots the output of pstats.ncl's average function by anomalies for
;    hours before/after all surges (i.e., 24 hr before, hour of, 24 hours after, etc.)
;
;    Like avgplot, this program also requires input of hoursBefore as it was in
;    the timegrab.ncl program as a positive number.  For simplicity, the 
;    anomalies are from the July/August All Hour mean, and as the input 
;    is surges from both July and August with no time specified, 
;    it is the most accurate mean that can be chosen.
;
;    The appropriate area is from 5 N to 47 N and 70 W to 130 W.
;
;    This version plots moisture flux vectors over color-filled theta-e.


; Load NCL graphics scripts and contributed functions
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/contributed.ncl"

;******************************************************************;
; SET hoursBefore TO BE THE SAME AS TIMEGRAB.NCL! (Positive value) ;
;******************************************************************;

hoursBefore=48

;*******************************************;
; Set file type, level, and variable number ;
;*******************************************;

filemid = "thetaE"
filemidu = "RS.clm.242"
filemidv = "RS.clm.243"
filemidout = "moistflx"


begin

; Use wc to count the number of files to be read in
numfilestemp = systemfunc("ls cathouravg."+filemid+".*.nc | wc -l")
numfiles=stringtointeger(numfilestemp)
fileArray = new(numfiles, string)
fileArrayu = new(numfiles, string)
fileArrayv = new(numfiles, string)
climFiles = new(numfiles, string)
climFilesu = new(numfiles, string)
climFilesv = new(numfiles, string)
indata = new((/277,349/),float)
indatau = new((/277,349/),float)
indatav = new((/277,349/),float)
cindata = new((/277,349/),float)
cindatau = new((/277,349/),float)
cindatav = new((/277,349/),float)
data = new((/277,349/),float)
datau = new((/277,349/),float)
datav = new((/277,349/),float)
plotc = new(6,graphic)
plotv = new(6,graphic)
plotbvect = new(6,graphic)
plotbcont = new(6,graphic)

hoursBefore = -hoursBefore

; Create the input file names
do i=0, numfiles-1
  fileArray(i) = "cathouravg."+filemid+"."+i+".nc"
  fileArrayu(i) = "cathouravg."+filemidu+"."+i+".nc"
  fileArrayv(i) = "cathouravg."+filemidv+"."+i+".nc"
end do

cf = "../../climodata/averages/JAAllAvg."+filemid+".nc"
climofile = addfile(cf, "r")
cf = "../../climodata/averages/JAAllAvg."+filemidu+".nc"
climofileu = addfile(cf, "r")
cf = "../../climodata/averages/JAAllAvg."+filemidv+".nc"
climofilev = addfile(cf, "r")

; Open a workspace to begin creating the plot
;wsplota = gsn_open_wks("ps", "./graphs/anomavgNARRa."+filemidout+".6panel")
;wsplotb = gsn_open_wks("ps", "./graphs/anomavgNARRb."+filemidout+".6panel")
wsplot = gsn_open_wks("ps", "./graphs/anomavgNARR."+filemidout+".6panel")

do i=0, numfiles-1, 8
  ; Tell user how far processing is
  print("Plotting hour "+(i+1)+"/"+numfiles)
  
  ; define data variables for within loop
  plotthis = new((/277,349/),float)
  plotthisu = new((/277,349/),float)
  plotthisv = new((/277,349/),float)
  vectorcolor = new((/277,349/),float)
  
  ; add input files for reading
  infile = addfile(fileArray(i), "r")
  infileu = addfile(fileArrayu(i), "r")
  infilev = addfile(fileArrayv(i), "r")
  
  ; Read in data from input file(s)
  indata = infile->avgNARR
  indatau = infileu->avgNARR
  indatav = infilev->avgNARR
  cindata = climofile->NARRfield
  cindatau = climofileu->NARRfield
  cindatav = climofilev->NARRfield
  
  ; calculate anomalies from the climatology
  data = indata-cindata
  datau = indatau-cindatau
  datav = indatav-cindatav
  
  ; Read lat-lon data (variable attributes) to indata
  nlat = 349-1
  mlon = 277-1
  data@lat2d = infile->gridlat_221
  datau@lat2d = infile->gridlat_221
  datav@lat2d = infile->gridlat_221
  lat2d = infile->gridlat_221
  data@lon2d = infile->gridlon_221
  datau@lon2d = infile->gridlon_221
  datav@lon2d = infile->gridlon_221
  lon2d = infile->gridlon_221

  ; Assign name and units to variables as needed
  datau@long_name = "Moisture Flux"
  datau@units = indatau@units
  data@long_name = "ThetaE"
  data@units = "K"
  
  ; Open cathourttest.thetaE.*.nc to get significance
  ;print("cathourttest."+filemid+"."+i+".nc")
  ttstopened = addfile("cathourttest."+filemid+"."+i+".nc","r")
  ; read ttest results
  ttst = ttstopened->ttestNARR
  ; make all data that is not significant look like missing data
  plotthis = mask(data,(ttst.lt.(0.9)),False)
  copy_VarMeta(data,plotthis)
  delete(ttst)
  
  ; Open cathourttest.RS.clm.242.*.nc to get significance
  print("cathourttest."+filemidu+"."+i+".nc")
  ttstopened = addfile("cathourttest."+filemidu+"."+i+".nc","r")
  ; read ttest results
  ttstu = ttstopened->ttestNARR
  thisu = where((ttstu.lt.(0.9)),0,datau)
  copy_VarMeta(datau,thisu)
  copy_VarMeta(datau,plotthisu)
  
  ; Open cathourttest.RS.clm.243.*.nc to get significance
  ;print("cathourttest."+filemidv+"."+i+".nc")
  ttstopened = addfile("cathourttest."+filemidv+"."+i+".nc","r")
  ; read ttest results
  ttstv = ttstopened->ttestNARR
  thisv = where((ttstv.lt.(0.9)),0,datav)
  copy_VarMeta(datav,thisv)
  copy_VarMeta(datav,plotthisv)
  
  ;fill vectorcolor with 0 for u-sig, 1 for both sig, and 2 for v-sig
  vecsize = dimsizes(thisu)
  vectorcolor@_FillValue = -999
  do x=0,vecsize(0)-1
    do y=0,vecsize(1)-1
      if(ismissing(thisu(x,y)).or.ismissing(thisv(x,y))) then
        vectorcolor(x,y) = -999
      else
        if((thisu(x,y).ne.0).and.(thisv(x,y).eq.0)) then
          vectorcolor(x,y) = 0
        else
          if((thisu(x,y).eq.0).and.(thisv(x,y).ne.0)) then
            vectorcolor(x,y) = 2
          else
            if((thisu(x,y).ne.0).and.(thisv(x,y).ne.0)) then
              vectorcolor(x,y) = 1
            else
              if((thisu(x,y).eq.0).and.(thisv(x,y).eq.0)) then
                vectorcolor(x,y) = -999
              end if
            end if
          end if
        end if
      end if
    end do
  end do

  ; convert back to grid-relative
  angle2 = sin(50.0*0.017453)*(lon2d+107.0)*(0.017453)
  sinx = sin(-angle2)
  cosx = cos(-angle2)
  plotthisu = (cosx*thisu)+(sinx*thisv)
  plotthisv = (-sinx*thisu)+(cosx*thisv)

  ; Change options for plots (res=True)
  resv = True
  resc = True
  resbvect = True
  resbcont = True
    
  ; set vector resources
  resv@vcRefMagnitudeF = 100000.0
  resv@vcRefLengthF = 0.008
  resv@vcGlyphStyle = "LineArrow"
  resv@vcRefAnnoOrthogonalPosF = -1.0
  resv@vcRefAnnoParallelPosF = 1.0
  resv@vcLineArrowThicknessF = 1.0
  resv@vcLineArrowHeadMinSizeF = 0.005
  resv@vcLineArrowHeadMaxSizeF = 0.005
  resv@vcVectorDrawOrder = "PostDraw"
  resv@vcMapDirection = False
  resv@vcScalarMissingValColor = -1
  ; blank strings above image
  resv@gsnCenterString = ""
  resv@gsnRightString = ""
  resv@gsnLeftString = ""
  
  ; copy settings from resv to resbvect
  copy_VarAtts(resv,resbvect)
  
  resv@vfXCStride = 4
  resv@vfYCStride = 4
  resv@sfXCStride = 4
  resv@sfYCStride = 4
  resv@sfMissingValueV = -999
  resv@vcMonoLineArrowColor = False
  resv@vcLevelSelectionMode = "ExplicitLevels"
  resv@vcLevels = (/0.5,1.5/)
  resv@vcLevelCount = 3
  resv@gsnSpreadColors = True
  
  resbvect@vfXCStride = 8
  resbvect@vfYCStride = 8
  resbvect@sfXCStride = 8
  resbvect@sfYCStride = 8
  resbvect@vcMonoLineArrowColor = True
  resbvect@vcLineArrowColor = (/0.8,0.8,0.8/)
  resbvect@vcRefAnnoOn = False
  
  ; Color fill on, contours off, line labels off, manual level selection 
  ; mode, spacing every 1.0 units, min value 270, max value 320, use 
  ; entire color spectrum
  resc@cnFillOn = True
  resc@cnLinesOn = False
  resc@cnLineLabelsOn = False
  resc@cnLevelSelectionMode = "ManualLevels"
  resc@cnInfoLabelOn = False
  resc@lbLabelPosition = "Center"
  resc@lbOrientation = "Vertical"
  resc@lbLabelDirection = "Down"
  resc@gsnSpreadColors = True
  
  ; Create the label bars for the colored contours
  resc@lbLabelFontHeightF = 0.020
  resc@lbTitleString = ""
  resc@lbLabelFont = 21
  resc@lbLabelAutoStride = False
  resc@lbLabelStride = 10
  resc@lbOrientation = "Vertical"
  resc@pmLabelBarOrthogonalPosF = 0.02
  resc@pmLabelZone = 2
  resc@lbBoxLinesOn = True
  resc@lbLabelBarOn = False
  ; blank strings above image
  resc@gsnCenterString = ""
  resc@gsnRightString = ""
  resc@gsnLeftString = ""
  
  ; copy settings from resc to resbcont
  copy_VarAtts(resc,resbcont)
  
  ;*******************************************************************;
  ;                     set contour max, min, interval                ;
                                                                      ;
  resc@cnLevelSpacingF = 0.2
  resc@cnMinLevelValF = -5
  resc@cnMaxLevelValF = 5
  
  resbcont@cnLevelSpacingF = 0.2
  resbcont@cnMinLevelValF = -5
  resbcont@cnMaxLevelValF = 5
                                                                      ;
  ;                     set contour max, min, interval                ;
  ;*******************************************************************;


  ; Draw map boundaries: coastlines, national boundaries, and states/provinces
  resbcont@mpAreaMaskingOn = False
  resbcont@mpFillOn = False
  resbcont@mpDataBaseVersion="Ncarg4_1"
  resbcont@mpOutlineBoundarySets = "National"
  resbcont@mpOutlineSpecifiers = (/"Canada:Provinces", "United States:States", "Mexico:States"/)

  ; Always display tick marks
  resbcont@pmTickMarkDisplayMode = "Always"
  ; Do not add cyclic reference point, maximize plot to workstation size,
  ; orient paper horizontally, 0.5 inch paper margin
  resbcont@gsnAddCyclic = False
  resbcont@gsnMaximize = True
  resbcont@gsnPaperOrientation = "portrait"
  resbcont@gsnPaperMargin = 0.5

  ; Define map projection properties
	resbcont@mpProjection = "LambertConformal"
	resbcont@mpLambertMeridianF = 253
 	resbcont@mpLambertParallel1F = 50
 	resbcont@mpLambertParallel2F = 50
 	
 	; Define display boundaries
 	resbcont@mpLimitMode = "Corners"
  resbcont@mpLeftCornerLatF = lat2d(0,nlat/6)
  resbcont@mpLeftCornerLonF = lon2d(0,nlat/6)
  resbcont@mpRightCornerLatF = lat2d(3*mlon/5,3*nlat/4)
  resbcont@mpRightCornerLonF = lon2d(3*mlon/5,3*nlat/4)
  ; map perimeter on, map fill on, no map rotation
  resbcont@mpPerimOn = True
  resbcont@mpCenterRotF = 0
  ; 10deg lat-lon grid spacing, show lat-lon grid, make lat-lon lines
  ; dashed
  resbcont@mpGridSpacingF = 10.0
  resbcont@mpGridAndLimbOn = True
  resbcont@mpGridLineDashPattern = 16
  
  ; NCDOverlay is false
  resv@tfDoNDCOverlay = False
  resc@tfDoNDCOverlay = False
  resbvect@tfDoNDCOverlay = False
  resbcont@tfDoNDCOverlay = False
  
  resv@lbLabelBarOn = False
  resbvect@lbLabelBarOn = False
  resc@lbLabelBarOn = False
  resbcont@lbLabelBarOn = False
  if(mod(i+1,8).eq.5)
;    resbcont@lbLabelBarOn = True
  end if
  if(mod(i+1,8).eq.1)
;    resc@lbLabelBarOn = True
  end if
  
  ; set colors for each plot
  resbcont@gsnSpreadColorStart = 32
  resbcont@gsnSpreadColorEnd = 72
  resc@gsnSpreadColorEnd = 102
  resv@gsnSpreadColorStart = 103
  resv@gsnSpreadColorEnd = 105

  ; Create image title
  if((hoursBefore+3*i).lt.0) then
    resbcont@tiMainString= "surge"+(hoursBefore + 3*i)+" hours"
  end if
  if((hoursBefore+3*i).eq.0) then
    resbcont@tiMainString= "surge arrival"
  end if
  if((hoursBefore+3*i).gt.0) then
    resbcont@tiMainString= "surge+"+(hoursBefore + 3*i)+" hours"
  end if
  
  ; don't draw the plot or advance the frame until all six images are drawn
  resv@gsnDraw = False
  resv@gsnFrame = False
  resc@gsnDraw = False
  resc@gsnFrame = False
  resbcont@gsnDraw = False
  resbcont@gsnFrame = False
  resbvect@gsnDraw = False
  resbvect@gsnFrame = False

  ; Define the colorscale for the plot and reverse it (R-W-B) and create plot
  gsn_define_colormap(wsplot,"BlWhRe")
  gsn_reverse_colormap(wsplot)
  newc1 = NhlNewColor(wsplot,0.5,0.5,0.5)
  newc2 = NhlNewColor(wsplot,0,0,0)
  newc3 = NhlNewColor(wsplot,0.5,0.5,0.5)
  
  ; plot background contours
  plotbcont(i/8) = gsn_csm_contour_map(wsplot,data(:,:),resbcont)
  ; plot significant contours
  plotc(i/8) = gsn_csm_contour(wsplot,plotthis(:,:),resc)
  
  ; plot background vectors
  plotbvect(i/8) = gsn_csm_vector(wsplot, datau(:,:), datav(:,:), resbvect)
  ; plot significant vectors
  plotv(i/8) = gsn_csm_vector_scalar(wsplot, plotthisu(:,:), plotthisv(:,:), vectorcolor(:,:), resv)
  
  ; overlay vectors on top of contours for subplot i/8
  overlay(plotbcont(i/8),plotc(i/8))
  overlay(plotbcont(i/8),plotbvect(i/8))
  overlay(plotbcont(i/8),plotv(i/8))

  ; Delete lots of variables
  delete(plotthis)
  delete(plotthisu)
  delete(plotthisv)
  delete(vectorcolor)
end do
  
; set panel resources
resP = True
; Create the label bar for the colored contours
resP@lbLabelFontHeightF = 0.020
resP@lbTitleString = ""
resP@lbLabelFont = 21
resP@lbLabelAutoStride = False
resP@lbLabelStride = 10
resP@lbOrientation = "Horizontal"
resP@pmLabelBarOrthogonalPosF = 0.0
resP@pmLabelZone = 2
resP@lbBoxLinesOn = True
resP@gsnPanelLabelBar = True
resP@txString = "Significant Theta-E and Moisture Flux Anomalies"
resP@txFontHeightF = 0.02
gsn_panel(wsplot,plotbcont,(/3,2/),resP)

delete(resP)
delete(resv)
delete(resc)
delete(resbcont)
delete(resbvect)

end