;=====Jilin=====

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/WRF_contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/cd_string.ncl"

begin

filepath = "5KM_ACONC.nc"  ;(ACONC file should be linked as "*.nc")
	a = addfile(filepath,"r")
	b = addfile("5KM_gridcro2d.nc","r")
	c = addfile("wrfout_d03.nc","r")
folderpath = "average/"
outfilepath=(/"NO2/","PM10/","PM25/","SO2/","O3/","CO/"/)

;=====First get the variables we will need

    PBLH_1 = c->PBLH(:,:,:)

    u10_1 = wrf_user_getvar(c,"U10",-1)    ; u at 10 m, mass point
    v10_1 = wrf_user_getvar(c,"V10",-1)    ; v at 10 m, mass point

    lat = b->LAT(0,0,:,:)
        dims = dimsizes(lat)
        print(dims)
    lon = b->LON(0,0,:,:)
        nx = dims(1)
        ny = dims(0)
        minlat = b->LAT(0,0,0,0)                ;y
        maxlat = b->LAT(0,0,ny-1,0)             ;ny=34
        minlon = b->LON(0,0,0,0)                ;x
        maxlon = b->LON(0,0,0,nx-1)             ;nx=52

    U10 = u10_1(:,1:ny,1:nx)
    V10 = v10_1(:,1:ny,1:nx)
    PBLH = PBLH_1(:,1:ny,1:nx)

;=====Average

    printVarSummary(U10)
    printVarSummary(PBLH)

    WNDSPD = wind_speed(U10(:,:,:),V10(:,:,:))

    PBLH_AVG = dim_avg_n(PBLH,0)
    WNDSPD_AVG = dim_avg_n(WNDSPD,0)

    wrf_smooth_2d( PBLH_AVG, 3 )
    wrf_smooth_2d( WNDSPD_AVG, 3 )

    PBLH_AVG@units = "m"
    PBLH_AVG@description = "Average PBLH"

    PBLH_AVG!0="lat"
    PBLH_AVG!1="lon"
    PBLH_AVG&lat=lat(:,0)
    PBLH_AVG&lon=lon(0,:)
    PBLH_AVG&lon@units="degrees_east"
    PBLH_AVG&lat@units="degrees_north"

;=====Set some Basic Plot options

  res = True
  res@gsnDraw  = False
  res@gsnFrame = False

  res@cnFillOn    = True         ;turn on color fill
  res@cnLinesOn   = False        ;turn off contour lines
  res@gsnMaximize = True

  shp_filename = ("jilin_shp/jilin_shp.shp")

;=====Set colorbar

  res@gsnAddCyclic = False
  res@gsnSpreadColors = True   ;use full range of color map
  res@lbOrientation = "vertical"  ;orientation
  res@lbBoxMinorExtentF        = 0.15
  res@lbLabelFontHeightF       = 0.010
  res@pmLabelBarOrthogonalPosF = -0.02
  res@vpYF                     = 0.95
  res@gsnMaximize = True

  ;res@tmXTOn = False         ;toptick
  ;res@tmYROn = False

  ;res@tmXBMode             = "Explicit"
  ;res@tmXBValues           = fspan(75,135,7)
  ;res@tmXBLabels           = TIME

;=====mpres

  mpres                        = True
  mpres@vpYF                   = 0.92
  mpres@vpXF                   = 0.12
  mpres@vpWidthF = 0.75 ; width 
  mpres@vpHeightF = 0.75 ; height

  mpres@gsnDraw  = False
  mpres@gsnFrame = False
  ;mpres@gsnAddCyclic           = False              ; regional data
  ;!!!!! any plot of data that is on a native grid, must use the "corners"
  ;method of zooming in on map.
  mpres@mpLimitMode       = "Corners"            ; choose range of map

  mpres@mpLeftCornerLatF  = lat(0,0)
  ;mpres@mpLeftCornerLatF  = lat(60,46)
  mpres@mpLeftCornerLonF  = lon(0,0)
  ;mpres@mpLeftCornerLonF  = lon(60,46)
  mpres@mpRightCornerLatF = lat(ny-1,nx-1)
  ;mpres@mpRightCornerLatF = lat(ny-1,nx-4)
  mpres@mpRightCornerLonF = lon(ny-1,nx-1)
  ;mpres@mpRightCornerLonF = lon(ny-1,nx-4)

  ;mpres@mpLimitMode = "LatLon"
  ;mpres@mpMinLatF               = 38.
  ;mpres@mpMaxLatF               = 48.
  ;mpres@mpMinLonF               = 116.
  ;mpres@mpMaxLonF               = 132.
  mpres@mpProjection          = "LambertConformal"
  mpres@mpLambertParallel1F   = 30
  mpres@mpLambertParallel2F   = 50
  mpres@mpLambertMeridianF    = 119      ;GRIDDESC

  mpres@pmTickMarkDisplayMode   = "Always"
  mpres@tmXTOn=False
  mpres@tmYROn=False
  ;res@cnFillDrawOrder         = "PreDraw"
  ;res@mpOutlineDrawOrder      = "PostDraw"      ; draw continental outline last
  mpres@mpDataBaseVersion       = "MediumRes"
  mpres@mpDataSetName           = "Earth..4"
  ;mpres@mpAreaMaskingOn         = True
  ;mpres@mpLandFillColor         = "white"
  ;mpres@mpInlandWaterFillColor  = "white"
  ;mpres@mpOceanFillColor        = "white"
  mpres@mpOutlineBoundarySets   = "NoBoundaries"
  mpres@mpGridAndLimbOn        =False
  mpres@mpGridLineDashPattern  = 2
  mpres@mpGeophysicalLineThicknessF = 2.0
  mpres@mpGridLineThicknessF        = 2.0
  mpres@mpLimbLineThicknessF        = 2.0
  mpres@mpNationalLineThicknessF    = 2.0
  mpres@mpUSStateLineThicknessF     = 2.0
  mpres@mpProvincialLineThicknessF  = 3.0

;=====SHP for GIS?

  gsres = True
  gsres@gsLineThicknessF = 3.0
  ;gsres@gsLineColor  = "grey"

;=====Plot Average PBLH

  print("Plot Average PBLH")
  cmap = RGBtoCmap("rgb.txt")
  wks_type="png"
  wks_type@wkWidth = 1800
  wks_type@wkHeight = 1800

  wks = gsn_open_wks("png","annual.PBLH")
  gsn_define_colormap(wks, cmap)

  pres = res
  ;pres@tiMainString = NO2@long_name + "  " + date(time)
  pres@gsnLeftString   = "Average PBLH"
  ;pres@gsnCenterString = "Forcast: "+ date2(time)
  pres@gsnRightString  = "PBL Height (m)"
  pres@txFontHeightF   = 0.015
  pres@cnLevelSelectionMode = "ExplicitLevels"  ;set explicit contour levels
  pres@cnLevels        = (/200,215,230,240,250,260,270,285,300,310,320,330,345,355,370,380,390,400,410,420,430,445,455,470,480,490,500,510,520,530,545,555,570,580,590,600,610,620,630,640,660,670,680,690,700,730,760,790,830,870,910,940,970,1000/)
  pres@cnExplicitLabelBarLabelsOn= True
  pres@lbLabelStrings  = (/"200",""," ","","","","","","300","","","","","","","","","400","","","","","","","","","500","","","","","","","","","600","","","","","","","","","700","","","","","","","","","1000"/)
  pres@lbBoxLinesOn    =False
  ;pres@cnLevels = (/0,35,75,115,150,250/)
  ;pres@cnFillColors = colors

;=====modified

;=====Plotting?

  plot = gsn_csm_contour(wks,PBLH_AVG,pres)
  map = gsn_csm_map(wks,mpres)

  overlay(map,plot)

  ;ploy = gsn_add_shapefile_polylines(wks,plot,shp_filename,gsres)
  ploy = gsn_add_shapefile_polylines(wks,map,shp_filename,gsres)

  ;draw(plot)
  draw(map)
  frame(wks)

  delete(pres)

end