绘制ACONC分布图

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/csm/contributed.ncl"

wrffile=addfile("/data1/loucx/Build_WRF/DATA/WRF/wrfout_d03_2020-07-01_00:00:00","r")
lat2d=wrffile->XLAT(0,3:185,3:185)
lon2d=wrffile->XLONG(0,3:185,3:185)
dims = dimsizes(lat2d)

DataDir1 = "/data1/loucx/CMAQ1/data/cctm/"
flist1 = systemfunc("ls "+DataDir1+"ACONC.*_d03")
Sfile1 = addfiles(flist1, "r")
ListSetType (Sfile1, "join")
isop=Sfile1[:]->ISOP
Isop=dim_avg_n(isop(:,:,0,:,:),0)*1000  ;Dimensions and sizes:[24] x [183] x [183]
Isop!0="day"
Isop!1="row"
Isop!2="col"
;do k=0,23
k=0
var_0=Isop(0,:,:)

   hour=k+8
 if(hour .ge. 24) then
   hour:=hour-24
 end if

     ;;;;;;;;;添加观测站;;;;;;
; f_h = asciiread("sitelist.csv",-1,"string")
; delim="," 
; lon=tofloat(str_get_field(f_h,5,delim))
; lat=tofloat(str_get_field(f_h,6,delim))
; idnum=str_get_field(f_h,1,delim)


     ;;;;;;;;;;;;;;;;;;;;;;;;;

 wks  = gsn_open_wks ("x11","day-average-isop"+hour+":00")       
 gsn_define_colormap(wks,"WhBlGrYeRe")

 res                        = True               ; plot mods desired
 res@gsnDraw                =False
 res@gsnFrame               =False
 res@gsnMaximize            = True
 res@tfDoNDCOverlay         = True ;;;;很可能是这行的原因
 res@mpDataSetName          = "Earth..4"
 res@mpDataBaseVersion = "MediumRes" 
 res@mpDataResolution = "Medium"
 res@mpProjection           = "LambertConformal"
 res@mpLambertParallel1F    = 30
 res@mpLambertParallel2F    = 60
 res@mpLambertMeridianF     = 118
 ;res@trGridType = "TriangularMesh"

 res@mpLimitMode            = "Corners"
 res@mpLeftCornerLatF       = lat2d(0,0)
 res@mpLeftCornerLonF       = lon2d(0,0)
 res@mpRightCornerLatF      = lat2d(dims(0)-1,dims(1)-1)
 res@mpRightCornerLonF      = lon2d(dims(0)-1,dims(1)-1)

; res@mpGridAndLimbOn = True ;是否绘制格网经纬线
; res@mpGridSpacingF = 10.
; res@mpGridLineDashPattern = 16
; res@mpGridLineThicknessF = 1.5
; res@mpGridLineColor = "grey"
; res@mpGridAndLimbDrawOrder = "PostDraw"
; res @mpOutlineBoundarySets = "Geophysical"
 
 res@mpGeophysicalLineThicknessF = 0.6 ;设置用于绘制地球物理边界轮廓的线的厚度
 res@pmTickMarkDisplayMode = "Always"
 res@tmXTOn = False ;不要绘制顶部轴刻度线
 res@tmYROn = False ;不要绘制右边轴刻度线
 res@tmXBLabelFontHeightF = 0.02
 res@cnFillOn = True
 res@cnLinesOn = False
 res@cnLevelSelectionMode = "ExplicitLevels"
;res@cnLevels=(/0,0.5,1,1.5,2,2.5,3,3.5,4,4.5,5,5.5,6,6.5,7,7.5,8/)
 res@cnLevels=(/0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2,2.1,2.2,2.3,2.4,2.5,2.6,2.7,2.8,2.9,3,3.1,3.2,3.3,3.4,3.5,3.6,3.7,3.8,3.9,4,4.1,4.2,4.3,4.4,4.5,4.6,4.7,4.8,4.9,5,5.1,5.2,5.3,5.4,5.5,5.6,5.7,5.8,5.9,6,6.1,6.2,6.3,6.4,6.5,6.6,6.7,6.8,6.9,7,7.1,7.2,7.3,7.4,7.5,7.6,7.7,7.8,7.9,8,8.1/)

 res@cnLineLabelsOn         = False
 res@lbLabelBarOn = True  ;开启图形各自的色标
 res@lbTopMarginF = 0.15
 res@pmLabelBarHeightF = 0.09
 res@pmLabelBarOrthogonalPosF = 0.02  ;上下移动
 res@lbLabelFontHeightF  = 0.0225
 res@lbTitleFontHeightF = 0.025
 res@lbTitleOffsetF = 0.12
 res@lbBoxLinesOn = False
 res@lbLabelStride =0.2

 res@tiMainString           = "ISOP-average"+hour+":00 (ppb)"
 res@tiMainFontHeightF      = 0.020              ; smaller title
 res@tiMainOffsetYF         = -0.005

 plres = True
 plres@gsnMaximize       = True
 plres@amJust   = "TopLeft"
 plres@gsnFrame         = False
 plres@gsnPanelDebug = False
 plres@gsnPanelBoxes = False
 plres@gsnPanelFigureStringsPerimOn = False
 plres@gsnPanelFigureStringsFontHeightF = 0.013
 plres@gsnPanelFigureStringsBackgroundFillColor = "transparent"

 res@lbLabelStride = 2
 dumstr1 = unique_string("marker1")
 plot = gsn_csm_contour_map(wks,var_0,res)
 shp1="/data1/loucx/Build_WRF/cnmap/cnhimap.shp"  
 lnres1        = True     
 lnres1@gsLineColor      ="black"
 lnres1@gsLineThicknessF = 3         ; 2x thickness  ---2 to 3.5
 shp_plot0 = gsn_add_shapefile_polylines(wks,plot,shp1,lnres1)
 gsn_panel(wks,(/plot/),(/1,1/),plres) 
 frame(wks)
;end do