;--------------------------------------------------
; These files are loaded by default in NCL V6.2.0 and newer
; 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/wrf/WRFUserARW.ncl"
;----------------------------------------------------------------------
; Main code
;----------------------------------------------------------------------
  begin

;**********************************************************************
;    Read the files and extract data
;**********************************************************************
; January

  temp_jan = new ((/120,66,31,10/),"float")

  yr = 0

  do yy = 2021,2030
  
  day = 0 

  dir  = "/home/sbasu05/projects/def-sauchyn/sbasu05/wrf_edmonton_exp/wrfout_"+yy+"/"

  do dd = 1,31
  
  if ( dd .lt. 10) then 

  filename = "wrfout_d03_"+yy+"-01-0"+dd+"_00:00:00"
  a = addfiles(dir + filename + ".nc","r")

  end if 

  if ( dd .ge. 10) then 

  filename = "wrfout_d03_"+yy+"-01-"+dd+"_00:00:00"
  a = addfiles(dir + filename + ".nc","r")

  end if   
   
  lat2d = wrf_user_getvar(a,"lat",0)
  lon2d = wrf_user_getvar(a,"lon",0)

  minlat = min(lat2d)   ; for use later
  maxlat = max(lat2d)
  minlon = min(lon2d)
  maxlon = max(lon2d)
  
  temp_exp = wrf_user_getvar(a,"T2",-1) 

  temp_mean_daily = dim_avg_n_Wrap(temp_exp,0)

  temp_jan(:,:,day,yr) = temp_mean_daily(:,:)

  day = day + 1

  end do
  
  yr = yr + 1

  end do

  printVarSummary(temp_jan)

  temp_mean_jan = dim_avg_n_Wrap(dim_avg_n_Wrap(temp_jan,3),2)

  temp_mean_jan@lat2d = lat2d   ; for plotting
  temp_mean_jan@lon2d = lon2d
  
  printVarSummary(temp_mean_jan)
  printMinMax(temp_mean_jan,True)

  delete(temp_mean_daily)
  delete(temp_exp)

;February

  temp_feb = new ((/120,66,28,10/),"float")

  yr = 0

  do yy = 2021,2030
  
  day = 0 

  dir  = "/home/sbasu05/projects/def-sauchyn/sbasu05/wrf_edmonton_exp/wrfout_"+yy+"/"

  do dd = 1,28
  
  if ( dd .lt. 10) then 

  filename = "wrfout_d03_"+yy+"-02-0"+dd+"_00:00:00"
  a = addfiles(dir + filename + ".nc","r")

  end if 

  if ( dd .ge. 10) then 

  filename = "wrfout_d03_"+yy+"-02-"+dd+"_00:00:00"
  a = addfiles(dir + filename + ".nc","r")

  end if   
   
  lat2d = wrf_user_getvar(a,"lat",0)
  lon2d = wrf_user_getvar(a,"lon",0)

  minlat = min(lat2d)   ; for use later
  maxlat = max(lat2d)
  minlon = min(lon2d)
  maxlon = max(lon2d)
  
  temp_exp = wrf_user_getvar(a,"T2",-1) 

  temp_mean_daily = dim_avg_n_Wrap(temp_exp,0)

  temp_feb(:,:,day,yr) = temp_mean_daily(:,:)

  day = day + 1

  end do
  
  yr = yr + 1

  end do

  printVarSummary(temp_feb)

  temp_mean_feb = dim_avg_n_Wrap(dim_avg_n_Wrap(temp_feb,3),2)

  temp_mean_feb@lat2d = lat2d   ; for plotting
  temp_mean_feb@lon2d = lon2d
  
  printVarSummary(temp_mean_feb)
  printMinMax(temp_mean_feb,True)

  delete(temp_mean_daily)
  delete(temp_exp)


;************************************************
; plotting parameters 
;************************************************
  wks  = gsn_open_wks("x11", "t2m_spatial_shapefile")        ; send graphics to PNG file

  plot = new(2,graphic)

  res                = True
  res@gsnFrame       = False
  res@gsnDraw        = False
  res@gsnLeftString  = ""
  res@gsnRightString = ""

  tf_res                       = res

  tf_res@lbLabelBarOn          = False 
 
  tf_res@pmLabelBarOrthogonalPosF = -0.005
  
  tf_res@sfXArray              = lon2d
  tf_res@sfYArray              = lat2d
  
  tf_res@mpGridAndLimbOn        = False
  tf_res@pmTickMarkDisplayMode  = "Always"           ; turn on tickmarks
  tf_res@mpGridLineDashPattern  = 5                  ; lat/lon lines dashed
  tf_res@mpGeophysicalLineThicknessF = 2.5
  tf_res@mpMinLatF             = 52.5     
  tf_res@mpMaxLatF             = 54
  tf_res@mpMinLonF             = 246
  tf_res@mpMaxLonF             = 248
  tf_res@mpCenterLonF          = 247
  tf_res@mpGridLatSpacingF     = 0.1
  tf_res@mpGridLonSpacingF     = 0.1

  tf_res@tmXBLabelFontHeightF  = 0.01

  tf_res@cnFillOn              = True
  tf_res@cnFillMode            = "RasterFill"   
  tf_res@cnLinesOn             = False 
  tf_res@cnFillPalette         = "MPL_coolwarm"
  tf_res@cnLevelSelectionMode  = "ManualLevels"
  tf_res@cnMinLevelValF        = 257
  tf_res@cnMaxLevelValF        = 295
  tf_res@cnLevelSpacingF       = 0.1
  tf_res@cnLineLabelsOn        = False

  plot(0) = gsn_csm_contour_map(wks,temp_mean_jan,tf_res)

  plot(1) = gsn_csm_contour_map(wks,temp_mean_feb,tf_res)

;---Attach the shapefile polylines using files read off gadm.org/country.

  canada_shp_name = "gadm36_CAN_3.shp"

  lnres                  = True
  lnres@gsLineColor      = "gray25"
  lnres@gsLineThicknessF = 0.5   

  can_id = gsn_add_shapefile_polylines(wks,plot(0),canada_shp_name,lnres)

;************************************************
; create panel
;************************************************
  resP = True

  resP@gsnPanelFigureStrings        = (/"Jan","Feb"/) ; add strings to panel

  resP@gsnPanelLabelBar             = True 
  resP@lbOrientation                = "Vertical"
  resP@lbLabelStride                = 10
  resP@lbBoxLinesOn                 = False


  resP@amJust                       = "BottomLeft"

  resP@gsnPanelFigureStringsPerimOn = False
  resP@gsnPanelFigureStringsBackgroundFillColor = -1

  gsn_panel(wks,plot,(/1,2/),resP)             ; now draw as one plot

  end