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/wrf/WRFUserARW.ncl"
;Xiaoming Hu xhu@ou.edu  2015 Nov

begin
layer = 0
  fileObs = systemfunc("ls *LT.csv.txt")

 simname = "chem3.7.1_NARR2d_CONUS_NLCD_NEI2011_MEGAN"
iday =  27 
if (iday.lt.10) then
folder = "wrf"+simname+".2011080"+iday+"00"
else
folder = "wrf"+simname+".201108"+iday+"00"
end if


do idomain =  2,2  

  files = systemfunc ("ls wrfout_d0"+idomain+"_2011-0*:00")

 do ifile = 6 ,  dimsizes(files)-1  ; obs LT, sim UTC  4 hours difference
  file_target = files(ifile) 
  f     = addfile(file_target+".nc","r")
  fmap  = f 
  Times_char = f->Times
   Times_char(0,13) = Times_char(0,4) 
   Times_char(0,16) = Times_char(0,4) 

;  ua  = wrf_user_getvar(f,"ua",0)   ; u on mass points
;  va  = wrf_user_getvar(f,"va",0)   ; v on mass points 
  x     = f->o3(0,layer,:,:)
  x = (/x*1000./)
  v     = f->V10(layer,:,:) 
  u     = f->U10(layer,:,:)
  ifile_obs = ifile-6  ; EST and UTC 
  figurename = "wrfout_d0"+idomain+"_o3OverlayObsAQSinfo_"+ifile_obs
;  figurename = "test"
  wks = gsn_open_wks("png" ,figurename)           ; ps,pdf,x11,ncgm,eps
  gsn_define_colormap(wks,"BlAqGrYeOrReVi200") ; select color map
      gsres               = True
      gsres@gsMarkerIndex = 16          ; circle at first
      gsres@gsMarkerThicknessF = 1
      gsres@gsMarkerSizeF = 0.01
      gsres@gsMarkerColor      = "black"
  tres       =  True
  tres@txFontHeightF = 0.02


  res                       = True             ; plot mods desired
  res@lbLabelFontHeightF = 0.018
  res@tmXBLabelFontHeightF = 0.018
  res@tmYLLabelFontHeightF = 0.018

  res@gsnFrame  = False
  res@gsnMaximize = True
  res@gsnPaperOrientation   = "portrait"
  res@gsnSpreadColors       = True             ; use full range of colormap
  res@cnFillOn              = True             ; color plot desired
  res@cnLinesOn             = False            ; turn off contour lines
  res@cnLineLabelsOn        = False            ; turn off contour labels
  res@cnLevelSelectionMode  = "ManualLevels"
  res@cnMaxLevelValF        = 90 
  res@cnMinLevelValF       =  10;1940  
  res@cnLevelSpacingF   =  5 

  res@lbLabelStride         =4 
;  res@lbLabelAutoStride = True ;
  res@tmYLLabelStride = 2
  res@tmXBLabelStride = 2

  WRF_map_c(fmap, res, 0)                   ; reads info from file
  res@mpOutlineBoundarySets = "AllBoundaries"
   res@vcRefAnnoOn = False
   res@vcRefMagnitudeF          = 8.             ; define vector ref mag
   res@vcRefLengthF             = 0.025          ; define length of vec ref
   res@vcRefAnnoOrthogonalPosF  = -1.            ; move ref vector
   res@vcRefAnnoParallelPosF  = 0.95            ; move ref vector
   res@vcMinDistanceF           = 0.025            ; larger means sparser
   res@vcLineArrowHeadMaxSizeF  = 0.0075          ; default: 0.05 (LineArrow), 0.012 (CurlyVector) 
   res@vcGlyphStyle             = "CurlyVector"   ; default: "LineArrow"  
   res@gsnScalarContour      = True               ; contours desired    

  res@tfDoNDCOverlay       = True     

  res@pmTickMarkDisplayMode = "Always"         ; turn on tickmarks
 res@tmXTOn = False            ; turn off top   labels 
 res@tmYROn = False            ; turn off right labels
; res@tiMainString     = "Valid " + chartostring(fmap->Times) + "(UTC) (Init: 2010080600; "+ifile+"hr forecast)"  
; res@gsnLeftString     = "T2 " + chartostring(fmap->Times) ;+ "(UTC) (Init: 2010080600; "+ifile+"hr forecast)"  
; res@gsnRightString    = "g kg~S~-1" 
 res@gsnLeftString = "O~B~3~N~"+ " "+chartostring(fmap->Times)

  res@cnInfoLabelOn       = True

  res@lbTitleOn        =  True                ; turn on title
;  res@lbTitleString    = "   ~C~  ~C~~V10~K"   ; "O~B~3, ~N~ppbv"             ; title string
  res@lbTitleString    = "  ~C~ ~C~  ~C~~V10~ppbv"   ; "O~B~3, ~N~ppbv"             ; title string

  res@lbTitlePosition  = "Right"              ; title position
  res@lbTitleFontHeightF= .018                ; make title smaller
  res@lbTitleDirection = "Across"             ; title direction
;  res@lbTopMarginF = -0.3
  res@lbTitleOffsetF = -0.01
  res@lbTitleJust = "CenterCenter"


  res@cnInfoLabelOrthogonalPosF = -0.04
  res@cnInfoLabelString   = "Min= $ZMN$ Max= $ZMX$"

  res@lbLabelFontHeightF = 0.018
  matrix = readAsciiTable(fileObs(ifile_obs),4,"float",0)
  string_xhu = fileObs(ifile_obs)
  string_xhu_char = stringtochar(string_xhu)
  string_xhu_show = charactertostring (string_xhu_char(5:19))
  res@gsnRightString    =  string_xhu_show 
  plot                 = gsn_csm_vector_scalar_map(wks,u,v,x,res)
  gsn_polymarker(wks,plot,-76.87647,39.056406,gsres)
;  gsn_text(wks,plot,"Research site",-76.87647,39.056406+0.2,tres)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  lat  = matrix(:,2)
  lon  = matrix(:,3)
  HourlyO3= matrix(:,1)

  indexbad = ind(HourlyO3.le.0)
   if (.not.all(ismissing(indexbad))) then
      HourlyO3(indexbad) = HourlyO3@_FillValue
   end if
   delete(indexbad);  = ind(HourlyO3.le.0)

  npts = dimsizes(lat)
  getvalues plot@contour
    "cnLevels" : levels
    "cnFillColors" : colors
  end getvalues
  num_distinct_markers = dimsizes(levels)+1        ; number of distinct markers
  dims=dimsizes(HourlyO3) 
  lat_new = new((/num_distinct_markers,dims/),float,-999)
  lon_new = new((/num_distinct_markers,dims/),float,-999)
  do j = 0, num_distinct_markers-1
    if (j.eq.0) then
      indexes = ind(HourlyO3.lt.levels(0))
    end if
    if (j.eq.num_distinct_markers-1) then
      indexes = ind(HourlyO3.ge.max(levels))
     print("color "+j+" indexes="+indexes)
    end if
    if (j.gt.0.and.j.lt.num_distinct_markers-1) then
      indexes = ind(HourlyO3.ge.levels(j-1).and.HourlyO3.lt.levels(j))
    end if

    if (.not.any(ismissing(indexes))) then
      npts_range = dimsizes(indexes)   ; # of points in this range.
      lat_new(j,0:npts_range-1) = lat(indexes)
      lon_new(j,0:npts_range-1) = lon(indexes)
     print("color "+j+"npts_range="+npts_range)
    end if
    delete(indexes)            ; Necessary b/c "indexes" may be a different
                               ; size next time.
  end do ; j = 0, num_distinct_markers-1
  gsres               = True
  gsres@gsMarkerIndex = 4          ; circle at first

  do j = 0, num_distinct_markers-1
    if (.not.ismissing(lat_new(j,0)))
      gsres@gsMarkerThicknessF = 1
      gsres@gsMarkerSizeF = 0.011
      gsres@gsMarkerColor      = "black"
      gsn_polymarker(wks,plot,lon_new(j,:),lat_new(j,:),gsres)
       gsres@gsMarkerIndex = 16          ; Use filled dots for markers.
      gsres@gsMarkerSizeF = 0.009
      gsres@gsMarkerColor      = colors(j)
      gsn_polymarker(wks,plot,lon_new(j,:),lat_new(j,:),gsres)
    end if
  end do

 delete(lat_new)
 delete(lon_new)
 delete(HourlyO3)
  delete(matrix)
 delete(lat)
 delete(lon)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  tres@txFontHeightF = 0.03
 ; gsn_text(wks,plot,"a",-102.4,40.8,tres)
  frame(wks)
; system("cp convert.ncl convert_temp.ncl")
; system("ex +g/figurename/s//"+figurename+"/g +wq convert_temp.ncl")
; system("~/tools/ncl_ncarg-5.2.1.Linux_x86_64_gcc346/bin/ncl convert_temp.ncl")
; system("convert -trim "+figurename+".eps /nsftor/xhu/public_html/HurricaneImpactonO3/WRFV3.4.1/BouLac/wrfchem3.4.1fnl2d_2011DFW.2011092000/"+figurename+".png")

;  if (idomain.eq.1)then
;   system("convert -trim -density 300 -resize 1000x1000 "+figurename+".eps /nsftor/xhu/public_html/HurricaneImpactonO3/WRFV3.7/BouLac/"+folder+"/"+figurename+".png")
;  else
 ;  system("convert -trim "+figurename+".eps                                /nsftor/xhu/public_html/HurricaneImpactonO3/WRFV3.7/BouLac/"+folder+"/"+figurename+".png")
;  end if

; system("rm "+figurename+".eps")

 print("finish plotting "+figurename+".eps")
 end do ; ifile
 delete(u)
 delete(v)
 delete(x)
 delete(files)
end do ; idomain
end