max_domain

; - Plotting WRF data
; - Drawing a time series plot
;----------------------------------------------------------------------
; These files are loaded by default in NCL V6.2.0 and newer
load "NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"

load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"

begin
;---Open file; substitute your own WRF output file here
fname = systemfunc("csh -c 'ls wrfout_"+"d04_"+"*00:00'")+".nc"
;print(fname)
;exit
f = addfiles(fname,"r")

ListSetType (f, "cat") ; concatenate (default)
;print(f)
;exit

xlat = f[0]->XLAT
xlon = f[0]->XLONG

; Read character variable Times
; Convert to units of "hours since" for plotting purposes
times = f[:]->Times
Times = wrf_times_c(f[:]->Times, 0) ; convert to "hours since"
; print(Time)
Time = str_sub_str(wrf_user_getvar(f,"times",-1),""," ")
do itime = 0, dimsizes(Time)-1
Time(itime) = systemfunc("date '+%Y-%m-%d%H:00:00' -d '"+Time(itime)+" 8 hour'")
end do

pm10 = f[:]->PM10 ;
pm25 = f[:]->PM2_5_DRY ; inded_01h_pre = wrf_user_getvar(input,"PREC_ACC_C",-1)+wrf_user_getvar(input,"PREC_ACC_NC",-1)
d01h_pre = wrf_user_getvar(f,"PREC_ACC_C",-1)+wrf_user_getvar(f,"PREC_ACC_NC",-1)
alist_pre = dim_max_n(d01h_pre,(/1,2/))
pm10_temp = pm10(:,0,:,:)
alist_pm10 = dim_max_n(pm10_temp,(/1,2/))
pm25_temp = pm25(:,0,:,:)
alist_pm25 = dim_max_n(pm25_temp,(/1,2/))
d_01h_CCN1 = wrf_user_getvar(f,"CCN1",-1)
d_01h_CCN1_temp = d_01h_CCN1(:,0,:,:)
;d_01h_CCN1M = max(d_01h_CCN1_temp)
d_01h_CCN1M = dim_max_n(d_01h_CCN1_temp,(/1,2/))
d_01h_CCN2 = wrf_user_getvar(f,"CCN2",-1)
d_01h_CCN2_temp = d_01h_CCN2(:,0,:,:)
;d_01h_CCN2M = max(d_01h_CCN2_temp)
d_01h_CCN2M = dim_max_n(d_01h_CCN2_temp,(/1,2/))
d_01h_CCN3 = wrf_user_getvar(f,"CCN3",-1)
d_01h_CCN3_temp = d_01h_CCN3(:,0,:,:)
;d_01h_CCN3M = max(d_01h_CCN3_temp)
d_01h_CCN3M = dim_max_n(d_01h_CCN3_temp,(/1,2/))
d_01h_CCN4 = wrf_user_getvar(f,"CCN4",-1)
d_01h_CCN4_temp = d_01h_CCN4(:,0,:,:)
;d_01h_CCN4M = max(d_01h_CCN4_temp)
d_01h_CCN4M = dim_max_n(d_01h_CCN4_temp,(/1,2/))
d_01h_CCN5 = wrf_user_getvar(f,"CCN5",-1)
d_01h_CCN5_temp = d_01h_CCN5(:,0,:,:)
;d_01h_CCN5M = max(d_01h_CCN5_temp)
d_01h_CCN5M = dim_max_n(d_01h_CCN5_temp,(/1,2/))
d_01h_CCN6 = wrf_user_getvar(f,"CCN6",-1)
d_01h_CCN6_temp = d_01h_CCN6(:,0,:,:)
;d_01h_CCN6M = max(d_01h_CCN6_temp)
d_01h_CCN6M = dim_max_n(d_01h_CCN6_temp,(/1,2/))
printVarSummary(d_01h_CCN6M)

;printVarSummary(voc)
;exit
;printVarSummary(o3)
;exit

;---find sites

;print(deglon(0))
;print(deglat(0))
;exit

opt = True

;---write time series into a table
;printVarSummary(so2)
;print(so2(:,22,loc(1,0),loc(0,0)))
;printVarSummary(no2)
;printVarSummary(pm10)
;printVarSummary(pm25)
;printVarSummary(o3)
;print(pm10(:,0,loc(1,0),loc(0,0)))
;exit
;;-----------

;write_table("s3_pm10.txt","w",alist_pm10,format)

varname = (/"pre", "ccn1","ccn2","ccn3","ccn4","ccn5","ccn6", "pm10", "pm25"/)
varform = (/"%12.2f", "%12.2f", "%12.2f", "%12.2f","%12.2f","%12.2f","%12.2f", "%12.2f", "%12.2f"/)
varlist = [/alist_pre,d_01h_CCN1M,d_01h_CCN2M,d_01h_CCN3M,d_01h_CCN4M,d_01h_CCN5M,d_01h_CCN6M,alist_pm10,alist_pm25/]
head = "BJT "+str_upper(str_concat(str_insert(varname,"",-12)))
write_table("d04_"+"chem_100.txt","w",[/head/],"%s")
dims = 1
text = ndtooned(Time)
;print(text)
;print(varform(0))
;print(varlist[0])

do ivar = 0, dimsizes(varname)-1
text = text+ndtooned(sprintf(varform(ivar),varlist[ivar]))
end do
write_table("d04_"+"chem_100.txt","a",[/text/],"%s")
delete([/varname,varform,varlist,text/])

;---Create plots
wks1 = gsn_open_wks("eps","WRF_time_series_pre")
wks2 = gsn_open_wks("eps","WRF_time_series_ccn3")
wks3 = gsn_open_wks("eps","WRF_time_series_ccn4")
wks4 = gsn_open_wks("eps","WRF_time_series_pm10")
wks5 = gsn_open_wks("eps","WRF_time_series_pm25")
plot1 = new(1,graphic)
plot2 = new(1,graphic)
plot3 = new(1,graphic)
plot4 = new(1,graphic)
plot5 = new(1,graphic)
gsn_define_colormap(wks1,"CBR_set3")
;colors = ispan(0,nsite-1,1)+2
;mkmode = ispan(0,nsite-1,1)+6
res = True ; plot mods desired
res@gsnDraw = False ; do not draw the plot (plots will be paneled)
res@gsnFrame = False
res@gsnMaximize = False ; maximize plot size
; Time@long_name = Time@units ; will be label for X-axis
res@tiMainString = "SO2"
res@vpWidthF = .8
res@vpHeightF = .5
res@xyMarkLineMode = "MarkLines"
res@xyLineThicknessF = 2.0 ; make thicker
res@tiYAxisString = "Precipitation (mm/h)"
res@trXMinF = 0
res@trXMaxF = 131
res@tmXBMode = "Explicit"
; res@tmXBValues = (/0,6,12,18,24,30,36,42,48,54,60,66,72,78,84,90,96,102,108,114,120/)
; res@tmXBLabels = (/"160504-00","06","12","18","0505-00","06","12","18","0506-00","06","12","18","0507-00","06","12","18","0508-00","06","12","18","0509-00","06"/)
res@tmXBValues = (/0,24,48,72,96/)
res@tmXBLabels = (/"Jun08","Jun09","Jun10","Jun11","Jun12"/)
res@tiXAxisString = ""

; res@vpHeightF = 0.4 ; change aspect ratio of plot
; res@vpWidthF = 0.70
;res@xyMarkers = mkmode
;res@xyMarkerColors = colors
res@xyMonoDashPattern = True ; all solid lines
;res@xyLineColors = (/"foreground","blue","red","brown"/)
;res@xyLineColors = colors

lgres = True
lgres@vpWidthF = .4
lgres@vpHeightF = .04
lgres@lgPerimOn = False
lgres@lgPerimColor = "black"
lgres@lgPerimThicknessF = 2.
lgres@lgItemType = "MarkLines"
lgres@lgMarkerSizeF = .01
lgres@lgLineThicknessF = 2.
lgres@lgLabelJust = "CenterCenter"
lgres@lgLabelFontHeightF = .05
lgres@lgLabelOffsetF = .1
amres = True
amres@amParallelPosF = -.4
amres@amOrthogonalPosF = .7

resP = True
resP@gsnPanelMainString = ""
resP@gsnPaperOrientation = "portrait"
resP@gsnMaximize = True
resP@gsnPanelLabelBar = False ; add common colorbar
resP@cnInfoLabelOn = False ; turn off cn info label
resP@gsnPanelRowSpec = True
resP@gsnPanelCenter = True
printVarSummary(Times)
printVarSummary(alist_pre)
plot1 = gsn_csm_xy(wks1,Times,(/alist_pre/),res)
resP= True
gsn_panel(wks1,plot1,(/1,1/),resP)

gsn_define_colormap(wks2,"CBR_set3")
;res@xyExplicitLegendLabels = site ; explicit labels created above

res@tiMainString = "CCN"
res@tiYAxisString = "CCN (/cm3)"
plot2 = gsn_csm_xy(wks2,Times,d_01h_CCN3M,res)
resP= True
gsn_panel(wks2,plot2,(/1,1/),resP)

gsn_define_colormap(wks3,"CBR_set3")
res@tiMainString = "CCN"
res@tiYAxisString = "CCN(/cm3)"
plot3 = gsn_csm_xy(wks2,Times,(/d_01h_CCN1M/),res)
resP= True
gsn_panel(wks3,plot3,(/1,1/),resP)

gsn_define_colormap(wks4,"CBR_set3")
res@tiMainString = "PM10"
res@tiYAxisString = "PM10 (ug/m3)"
;res@tiXAxisString = "Date"
plot4 = gsn_csm_xy(wks4,Times,(/alist_pm10/),res)
resP= True
gsn_panel(wks4,plot4,(/1,1/),resP)

gsn_define_colormap(wks5,"CBR_set3")
res@tiMainString = "PM2.5"
res@tiYAxisString = "PM2.5 (ug/m3)"
plot5 = gsn_csm_xy(wks5,Times,(/alist_pm10/),res)
resP= True
gsn_panel(wks5,plot5,(/1,1/),resP)
end