MDA8 有/无生物ISOP的O3差值

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"
VALUE=new((/31,17,174,168/),float)
LUE=new((/31,17,174,168/),float) 
VALUE2=new((/174,168/),float)
VALUE3=new((/174,168/),float)
MAX1=new((/31/),float)
MAX2=new((/31/),float)
wrffile=addfile("/data1/loucx/COPY/WRF/WRF/run/wrfout_d03_2020-08-01_00:00:00","r")
lat2d=wrffile->XLAT(0,3:176,3:170)
lon2d=wrffile->XLONG(0,3:176,3:170)
dims = dimsizes(lat2d)
DataDir1 = "/data1/loucx/COPY/CNAQ/data1/cctm/"
DataDir2 = "/data1/loucx/COPY/CNAQ/data/cctm/"
flist1 = systemfunc("ls "+DataDir1+"ACONC.*_d03")
flist2 = systemfunc("ls "+DataDir2+"ACONC.*_d03")
Sfile1 = addfiles(flist1, "r")
Sfile2 = addfiles(flist2, "r")
ListSetType (Sfile1, "cat")
ListSetType (Sfile2, "cat")
o31=Sfile1[:]->O3 ;[TSTEP | 744] x [LAY | 1] x [ROW | 174] x [COL |168]
O31=o31(:,0,:,:) ;[TSTEP | 744] x [ROW | 174] x [COL | 168]
do r=0,173
 do c=0,167
     ;7.1
     k=7  ;位置
    do i=0,9  ;ID
     j=i+7
     VALUE(0,k,r,c)=avg(O31(i:j,r,c))
     k=k+1
    end do
  end do
end do     

;7.2-7.30
do r=0,173
 do c=0,167
     da=1
    do i=17,689,24 ;day
      k=0
      do x=i,i+16   
        j=x+7
        VALUE(da,k,r,c)=avg(O31(i:j,r,c))
        k=k+1
      end do
      da=da+1
    end do
  end do
end do
;7.31
do r=0,173
 do c=0,167
     k=0  ;位置
    do i=713,728  ;ID
     j=i+7
     VALUE(30,k,r,c)=avg(O31(i:j,r,c))
     k=k+1
    end do
  end do
end do   
;MDA8
do r=0,173
 do c=0,167
  do t=0,30
   MAX1(t)=max(VALUE(t,:,r,c))
  end do
  VALUE2(r,c)=avg(MAX1(:))*48000/22.4  ;[174] x [168] *48000/22.4 ug/m-3
 end do
end do

o32=Sfile2[:]->O3 ;[TSTEP | 744] x [LAY | 1] x [ROW | 174] x [COL |168]
O32=o32(:,0,:,:) ;[TSTEP | 744] x [ROW | 174] x [COL | 168]
do r=0,173
 do c=0,167
     ;7.1
     k=7  ;位置
    do i=0,9  ;ID
     j=i+7
     LUE(0,k,r,c)=avg(O32(i:j,r,c))
     k=k+1
    end do
  end do
end do     

;7.2-7.30
do r=0,173
 do c=0,167
     da=1
    do i=17,689,24 ;day
      k=0
      do x=i,i+16   
        j=x+7
        LUE(da,k,r,c)=avg(O32(i:j,r,c))
        k=k+1
      end do
      da=da+1
    end do
  end do
end do
;7.31
do r=0,173
 do c=0,167
     k=0  ;位置
    do i=713,728  ;ID
     j=i+7
     LUE(30,k,r,c)=avg(O32(i:j,r,c))
     k=k+1
    end do
  end do
end do   
;MDA8
do r=0,173
 do c=0,167
  do t=0,30
   MAX2(t)=max(LUE(t,:,r,c))
  end do
  VALUE3(r,c)=avg(MAX2(:))*48000/22.4  ;[174] x [168] *48000/22.4 ug/m-3
 end do
end do

var_0=VALUE2(:,:)-VALUE3(:,:)
wks  = gsn_open_wks ("x11","Avg8-MDA8")       
 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@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,2,4,6,8,10,12,14,16,18,20,22/)

 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           = "Avg8-MDA8"
 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)