;*******************************************
; Written by Yuqiang Zhang, UNC, 2014
; Modifie by JHuang 20151013
; Concepts illustrated:
;   - Plotting CMAQ data; also works for data using LCC project, like WRF and other RCM
;   - Drawing filled contours over a Lambert Conformal map
;   - Drawing U.S. states
;********************************************************************
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/csm/contributed.ncl"
;load ".hluresfile"

begin
;**************************
;reading LAT and LON data
;**************************
   AP = (/"LAX","ORD","ATL"/)
   COL = (/21,96,108/)
   ROW = (/45,65,40/) 

l = 1
   LON_1 = COL(l)-5
   LON_2 = COL(l)+5
   LAT_1 = ROW(l)-5
   LAT_2 = ROW(l)+5

   a = -0.004
   b = 0.004
   c = 0.0004
   ;read the 2D lon&lat
   lonlatfile = "/nas/depts/003/cep-emc/proj/FAA-D2025/joey/domain/Pradeepa_GRID/GRIDCRO2D_2004335_D2025_36km"
   flonlat = addfile(lonlatfile,"r")
   lon2d_temp = flonlat->LON
   lat2d_temp = flonlat->LAT
   printVarSummary(lon2d_temp)
   lon = lon2d_temp(0,0,LAT_1:LAT_2,LON_1:LON_2) 
   lat = lat2d_temp(0,0,LAT_1:LAT_2,LON_1:LON_2)
   lon@units = "degrees_east"
   lat@units = "degrees_north"
   printVarSummary(lon)
   lat2 = lat2d_temp(0,0,ROW(l),COL(l))
   lon2 = lon2d_temp(0,0,ROW(l),COL(l))
;**************************
;reading CMAQ data
;**************************

; reading scenario 1 (S1) data
CMAQ = "v5.0.2"
case = (/"airc_base","airc_sens_6","airc_bkgd"/)
   nameofpoll = (/"PMIJ","AECIJ","ASO4IJ","AOCIJ","ANH4IJ","ANO3IJ"/)
   name_title = (/"PM~B~2.5~N~","EC","SO~B~4~N~","OC","NH~B~4~N~","NO~B~3~N~"/)
   npoll = dimsizes (nameofpoll)
i = 0
   wks = gsn_open_wks("png","delta_"+AP(l)+"_"+nameofpoll(i)+"_annual_2005")          ; open a workstation
   gsn_define_colormap(wks, "BlueWhiteOrangeRed") ;"BlAqGrYeOrReVi200");(wks,"gui_default")  "BlAqGrYeOrRe" ; choose colormap
   plots = new(1, graphic)
k = 0
   file_data_base = "/nas/depts/003/cep-emc/data/FAA-D2025/ESF_test/ESF_"+CMAQ+"/"+case(2)+"_2005_post/"+case(2)+".ACONC.2005.combine.yave"  
   file_data_case = "/nas/depts/003/cep-emc/data/FAA-D2025/ESF_test/ESF_"+CMAQ+"/"+case(k)+"_2005_post/"+case(k)+".ACONC.2005.combine.yave"  
   
   print(file_data_base)
   print(file_data_case)

   f_1W = addfile(file_data_base,"r")
   f_2W = addfile(file_data_case,"r")

   POLLN=nameofpoll(i)
   POLLN_base = f_1W->$POLLN$
   POLLN_case = f_2W->$POLLN$
   att =  getvaratts(POLLN_base)
   unit = POLLN_base@$att(1)$

   POLL_diff_temp = (POLLN_case - POLLN_base)
   printVarSummary(POLL_diff_temp)

;domain average into two dimensions
   POLL_diff_1 = dim_avg_n_Wrap(POLL_diff_temp,(/0,1/))
   POLL_diff = POLL_diff_1(LAT_1:LAT_2,LON_1:LON_2)
   printVarSummary(POLL_diff)
   MA = max(POLL_diff)
   MI = min(POLL_diff)
   AVE = avg(POLL_diff)

;assign lon and lat attribution to data   
   POLL_diff@lon2d=lon
   POLL_diff@lat2d=lat


;get the dimensions
   dims=dimsizes(POLL_diff)                            ; get dimensions
   print(dims)
   ndims = dimsizes(dims)
   nlon=dims(ndims-1)                                ; assign # lat/lon points
   nlat=dims(ndims-2)
   print(nlon)
   print(nlat)

;**************************
;plotting
;**************************
   
; to define the parameters that needed for the plotting
   res                 = True                    ; plot mods desired
   res@cnLinesOn       = False                   ; default is true;False, no contour lines will appear in the contour plot
   res@cnLineLabelsOn  = False                   ; turn off contour lines
   res@cnFillOn        = True                    ; color plot desired
;   res@cnFillPalette   = ""                ; choose color map
   res@cnLevelSelectionMode = "ManualLevels"    ; manually set the contour levels with the following 3 resources
   res@cnMinLevelValF  = a                    ; set the minimum contour level
   res@cnMaxLevelValF  = b                 ; set the maximum contour level
   res@cnLevelSpacingF = c                    ; set the interval between contours
;   res@cnLevelSelectionMode = "ExplicitLevels" ; set explicit contour levels 
;   res@cnLevels = log10((/-1e8,-1e6,-1e4,-1e2,0,1e2,1e4,1e6,1e8 /)) ; set levels 

; !!!!! any plot of data that is on a native grid, must use the "corners"
; method of zooming in on map.   
   res@mpLimitMode       = "Corners"            ; choose range of map
   res@mpLeftCornerLatF  = lat(0,0)
   res@mpLeftCornerLonF  = lon(0,0)
   res@mpRightCornerLatF = lat(nlat-1,nlon-1)
   res@mpRightCornerLonF = lon(nlat-1,nlon-1)

; The following 4 pieces of information are REQUIRED to properly display
; data on a native lambert conformal grid.  This data should be specified
; somewhere in the model itself.      
   res@mpProjection        = "LambertConformal"
   res@mpLambertParallel1F = 33
   res@mpLambertParallel2F = 45
   res@mpLambertMeridianF  = -97


; usually, when data is placed onto a map, it is TRANSFORMED to the specified
; projection. Since this model is already on a native lambert conformal grid,
; we want to turn OFF the tranformation.
   res@tfDoNDCOverlay = True
   res@pmTickMarkDisplayMode = "Always"                 ; turn on tickmarks
   res@tmXTOn = False            ; turn off top   labels 
   res@tmYROn = False             ; turn off right labels
   res@tmXBOn = True            ; turn off bottom lables
   res@tmYLOn = True            ; turn off left lables

   res@mpGeophysicalLineColor = "Black"                   ; color of continental outlines
   res@mpPerimOn              = True                    ; draw box around map
   res@mpGridLineDashPattern  = 2                       ; lat/lon lines as dashed
   res@mpGridLineColor        = "Black"                 ; GridLine color
   res@mpOutlineBoundarySets = "GeophysicalAndUSStates" ; add state boundaries
   res@mpUSStateLineColor    = "Black"


; panel plots
  res@gsnDraw                = False     ; do not draw
  res@gsnFrame               = False     ; do not advance frame
  res@lbLabelBarOn = True              ;trun on labels
  res@lbLabelFontHeightF = 0.02

;  res@lbLabelBarOn = False              ;delete all labels
;  res@cnInfoLableOn = False              ;delete all labels
; remove variable label
   res@gsnLeftString = ""
   res@gsnRightString = ""


; Define tickmarks for the first two plots
;for the colorbar or called "tickmarks" in NCL
   res@cnLevelSelectionMode  = "ManualLevels"
  
;create the map
; Working on titles
   res@tiMainString = name_title(i)+" "+AP(l)+"_"+case(k)+"- base ~F33~m~F21~g m~S~-3~N~"
;   res@tiMainString = "UFP_num_"+case(k)+"- base #/m~S~-3~N~" ;~F33~m~F21~g m~S~-3~N~"
   res@tiMainFontHeightF = 0.02
   res@gsnCenterString = "max:"+sprintf("%5.1e",MA)+", min:"+sprintf("%5.1e",MI)+", mean:"+sprintf("%5.1e",AVE)   

;plotting WS2-S1
   plots(k) = gsn_csm_contour_map(wks,POLL_diff,res)
   mkres               = True
   mkres@gsMarkerIndex = 16     ; Filled circle
   mkres@gsMarkerSizeF = 5
   print(lon2)
   print(lat2)
   gsn_add_polymarker(wks,plots(k),lon2,lat2,mkres)

; create panel
;=============================================
  pres                         = True
  pres@gsnPanelRowSpec         = True      ;tell panel what order to plot
  gsn_panel (wks, plots, (/1/), pres)
  ;pres@gsnPanelCenter          = True      ;put the plot in center
  ;gsn_panel ()
  print("-------------------------------------------------------------------------")
  print("diff_"+nameofpoll(i))
  print("-------------------------------------------------------------------------")
end