;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; corrmap_t2.ncl
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

begin

; Ensemble members
  column = (/"ck6m","ck6y","cktm","ckty","cn6m","cn6y","cntm","cnty","ct6m","ct6y","cttm","ctty",\
             "rk6m","rk6y","rktm","rkty","rn6m","rn6y","rntm","rnty","rt6m","rt6y","rttm","rtty"/)
  numColumn  = dimsizes(column)
; Regions we care about
  row = (/"PacificNW","PacificSW","Southwest","GreatBasin","NRockies",\
          "SRockies","Mezquital","NPlains","CPlains","SPlains","Prairie","DeepSouth","Southeast",\
          "GreatLakes","Appalachia","MidAtlantic","NorthAtlantic","EastBoreal"/)
; Regions as we want them to appear on the labels
  row_label = (/"Pacific NW","Pacific SW","Southwest","Great Basin","North Rockies",\
                "South Rockies","Mezquital","North Plains","Central Plains","South Plains","Prairie",\
                "Deep South","Southeast","Great Lakes","Appalachia","Mid Atlantic",\
                "North Atlantic","East Boreal"/)
  numRow = dimsizes(row)

; Define the directory where our data is
  dir = ""
; Read in the reanalysis data
  cfsr = addfile(dir+"t2_daily_Regional_Data_jja_cfsr.nc","r")
; Create a new variable for the T2 diff between each model and cfsr
  t2_diff = new((/numRow,numColumn/),float)
; Some do loops that loop through each model (column) and region (row)
  do i=0,numColumn-1
    print("Ensemble member: "+column(i))
    do j=0,numRow-1
      print("	Region: "+row_label(j))
      var = "T2MAX_"+row(j)  ; Name of the variable we want
      data = addfile(dir+"t2_daily_Regional_Data_jja_"+column(i)+".nc","r") ; Read in model data
      t2 = data->$var$  ; Model T2MAX
      t2_cfsr = cfsr->$var$-273.15  ; Reanalysis T2MAX
      t2_diff(j,i) = avg(t2)-avg(t2_cfsr)  ; Average both T2MAX's across all dims and take diff
      delete([/var,data,t2_cfsr,t2/])  ; Delete variables
    end do  ; End row do loop
  end do  ; End column do loop

; Because our rows and columns are defined by text we need to give them numbers for plotting
  XB = new (numColumn,float)
  XB(0) = -0.5
  do xx = 1,numColumn-1
    XB(xx) = XB(xx-1) + 1.
  end do
  YL = new (numRow,float)
  YL(0) = -0.5
  do yy = 1,numRow-1
    YL(yy) = YL(yy-1) + 1.
  end do

; Starting the plotting!
  wks = gsn_open_wks("X11","corrmap_t2max")
  gsn_define_colormap(wks,"BlueWhiteOrangeRed")

; Setting lots of resources for the plot
  res = True
  res@gsnMaximize = True
  res@gsnDraw = False
  res@gsnFrame = False
  res@gsnSpreadColors = True
  res@vpWidthF = 0.80      ;Sizing the plot so the squares are actually square
  res@vpHeightF = 0.60
  res@vpXF = 0.01
  res@vpYF =0.95

  res@cnFillOn = True
  res@cnLinesOn = False
  res@cnLineLabelsOn = False
  res@cnFillMode = "CellFill" ; A raster type fill
  res@cnCellFillEdgeColor = "black" ; Default edges are transparent
  res@cnLevelSelectionMode = "ManualLevels"
  res@cnLevelSpacingF = 0.5
  res@cnMinLevelValF = -4.0
  res@cnMaxLevelValF = 4.0

  res@tiMainString = "T2MAX Model-CFSR JJA"
  res@lbLabelAutoStride = True
  res@lbTitleDirection = "Across"
  res@lbTitleString = "~F34~0~F~C" ; Degree symbol
  res@lbTitlePosition = "Right"

  ; Lots of messing with tick marks and labels to make it look
  ;  exactly how I want
  res@tmXBMode = "Explicit"
  res@tmXBLabelAngleF = 90.
  res@tmXBLabelJust = "TopRight"
  res@tmXBValues = XB(:)
  res@tmXBLabels = column(:)
  res@tmXBLabelFontHeightF = 0.015
  res@tmYLMode = "Explicit"
  res@tmYLValues = YL(:)
  res@tmYLLabels = row_label(:)
  res@tmYLLabelJust = "BottomRight"
  res@tmYLLabelFontHeightF = 0.015
  res@tmXBMajorLengthF = 0.0
  res@tmXTMajorLengthF = 0.0
  res@tmYLMajorLengthF = 0.0
  res@tmYRMajorLengthF = 0.0

  ; Some magic had to be done here to make sure the raster plots
  ;  correctly. The edges get cut off so we need to add an additional
  ;  element to our coordinate arrays. They must be one element longer
  ;  than the data.
  res@sfXArray = fspan(-0.5,23,numColumn+1)
  res@sfYArray = fspan(-0.5,17,numRow+1)

; Make the plot!
  plot = gsn_csm_contour(wks,t2_diff,res)

; Draw the plot
  draw(plot)
  frame(wks)

end