load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"
load "~/res_wrf.ncl"
begin
wks = gsn_open_wks("pdf" ,"terra_l3_test")
gsn_define_colormap(wks, "BlAqGrYeOrRe")
plot = new(2,graphic)

dir1 = "/HOME/ustc_czh_4/WORKSPACE/Model/BASE/data/"
fil = addfile(dir1+"optical.nc","r")
lon = wrf_user_getvar(fil,"XLONG",0)
lat = wrf_user_getvar(fil,"XLAT",0) 
dim = dimsizes(lat) 
  tauaer2_ari = wrf_user_getvar(fil,"TAUAER2",-1) ; 400nm aod
  tauaer3_ari = wrf_user_getvar(fil,"TAUAER3",-1) ; 600nm aod
  ; Sum up the AOD at different height to get the total column AOD
  aod2_ari = dim_sum_n(tauaer2_ari(:,:,:,:),1)
  aod3_ari = dim_sum_n(tauaer3_ari(:,:,:,:),1)
  delete([/tauaer2_ari,tauaer3_ari/])
  ;Interpolate AOD to 550 nm
  ;Angstrom exponent to compute 550nm aod
  ae_ari = -1*(log(aod2_ari/aod3_ari)/log(0.4/0.6))
  var_ari = aod3_ari*((0.55/0.6)^ae_ari)
  delete([/aod2_ari,aod3_ari,ae_ari/])
  var_ari_0am = var_ari(0::8,:,:)
  printVarSummary(var_ari_0am)
  var_ari_3am = var_ari(1::8,:,:)
  var_ari_6am = var_ari(2::8,:,:)
  var_terra = (var_ari_3am + var_ari_0am)/2
  var_terra@_FillValue = default_fillvalue("float")
  printVarSummary(var_terra)
; --- WRF boundaries
  minLatWrf  = min(lat)    
  maxLatWrf  = max(lat)
  minLonWrf  = min(lon)
  maxLonWrf  = max(lon)
dir2 = "/HOME/ustc_czh_4/WORKSPACE/Analysis/observation/modis/terra/"
fils = systemfunc("cd "+ dir2 + "; ls MOD08_D3.A2015*")
nf = dimsizes(fils)
aod_terra_wrf = new((/nf,dim(0),dim(1)/),"float")
do ifl = 0, nf-1
  fil0 = addfile(dir2+fils(ifl),"r")
  aod_terra = short2flt(fil0->AOD_550_Dark_Target_Deep_Blue_Combined_Mean)
  aod_terra@_FillValue = default_fillvalue("float")
  aod_terra&XDim_mod08 = fspan(-179.5,179.5,360)
  aod_terra&YDim_mod08 = fspan(89.5,-89.5,180)
  aod_terra&XDim_mod08@units = "degrees_east"
  aod_terra&YDim_mod08@units = "degrees_north"
  ;aod_terra_reg = aod_terra({minLatWrf:maxLatWrf},{minLonWrf:maxLonWrf})
; --- Regrid
  ;aod_terra_wrf(ifl,:,:) = rgrid2rcm_Wrap(aod_terra_reg&YDim_mod08, aod_terra_reg&XDim_mod08,\
  ;                aod_terra_reg, lat, lon,1) 
  aod_terra_wrf(ifl,:,:) = rgrid2rcm_Wrap(aod_terra&YDim_mod08, aod_terra&XDim_mod08,\
                   aod_terra, lat, lon,1) 
  printVarSummary(aod_terra_wrf)
  printMinMax(aod_terra_wrf,True)
  var_terra(ifl,:,:) =where(ismissing(aod_terra_wrf(ifl,:,:)),default_fillvalue("float"),var_terra(ifl,:,:))
  
  N1 = num(ismissing(aod_terra_wrf(ifl,:,:)))
  N2 = num(ismissing(var_terra(ifl,:,:)))
  print("obs missing value number is "+N1+"; model missing value number is "+N2)
end do

  aod_terra_wrf_ave = dim_avg_n(aod_terra_wrf,0)
  var_terra_ave = dim_avg_n(var_terra,0)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;        plot
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  res@lbBoxLinesOn          = False
  res@lbLabelBarOn = False
  res@cnLevels    = (/0.1,0.2,0.3,0.4,0.6,0.8,1.0,1.2/) ; set levels
  res@cnFillColors = (/2,12,17,34,60,73,81,89,101/)     ; set the colors to be used
  pltres@PlotTitle =  "(a) Terra"
  contour                   = wrf_contour(fil,wks,aod_terra_wrf_ave,res)
  plot(0) = wrf_map_overlays(fil,wks,contour,pltres,mpres)
  pltres@PlotTitle =  "(b) WRF-Chem_terra"
  contour                   = wrf_contour(fil,wks,var_terra_ave,res)
  plot(1) = wrf_map_overlays(fil,wks,contour,pltres,mpres)
  resP@gsnPanelLabelBar  = True
  resP@lbBoxLinesOn          = False
  resP@txString             = "AOD at 550nm"
  gsn_panel(wks,plot,(/2,1/),resP)

end