; Xiaoming Hu Sept. 19, 2011
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"

undef("add_map_tickmarks")
function add_map_tickmarks(wks,plot,res)
local res2, bres, vpx, vpy, vpw, vph, xndc, yndc, npts, n, j, nlat, \
nlon, delta, bot_lon, top_lon, lft_lat, rgt_lat, xblabels, xbvalues, \
xtlabels, xtvalues, yllabels, ylvalues, yrlabels, yrvalues, xfix, \
xlat, xlon, yfix, annoid, anno_str
begin
;---Make a copy of the original resource list.
  res2 = res

;---Retrieve edges of plot in NDC space.
  getvalues plot
     "vpXF"      :  vpx
     "vpYF"      :  vpy
     "vpWidthF"  :  vpw
     "vpHeightF" :  vph
  end getvalues

;---Turn off tickmarks associated with map. We want to add our own.
  setvalues plot
    "pmTickMarkDisplayMode" : "Never"
  end setvalues

;---Initialize resources for tickmark plot. User shouldn't change these.
  bres                          = True
  bres@vpXF                     = vpx
  bres@vpYF                     = vpy
  bres@vpWidthF                 = vpw
  bres@vpHeightF                = vph
  bres@trXMinF                  = vpx
  bres@trXMaxF                  = vpx + vpw
  bres@trYMinF                  = vpy - vph
  bres@trYMaxF                  = vpy
  bres@tmEqualizeXYSizes        = True

;---This resource the user can change in main code if desired.
  bres@gsnTickMarksPointOutward = get_res_value(res2,"gsnTickMarksPointOutward",True)

;
; NDC Points to scan on X and Y axes. These arrays will be used to
; find the closest NDC pair that gets us close to the location where
; we want a tickmark.
;
  npts = 100000   ; Increase to get closer match for tickmarks
  xndc = fspan(vpx,vpx+vpw,npts)
  yndc = fspan(vpy-vph,vpy,npts)

  n    = dimsizes(yndc) 
  xfix = new(n,float)
  yfix = new(n,float)
  xlon = new(n,float)
  xlat = new(n,float)
  delta = 0.001

;---Left axis tickmarks
  if(isatt(res2,"tmYLValues")) then
    lft_lat    = get_res_value(res2,"tmYLValues",-1)
    nlat       = dimsizes(lft_lat)
    ylvalues = new(nlat,float)
    yllabels = new(nlat,string)
    xfix  = vpx + 0.0001 ; Just a smidge into the plot to make sure we don't
                         ; get missing values returned.
;
; Loop across each left latitude value that we want a tickmark for,
; and try to find the closest X,Y NDC coordinate pair along this axis.
;
    NhlNDCToData(plot,xfix,yndc,xlon,xlat)
    do j=0,dimsizes(lft_lat)-1
      NhlNDCToData(plot,xfix,yndc,xlon,xlat)
      ii = minind(fabs(xlat-lft_lat(j)))
      if(.not.any(ismissing(ii)).and.fabs(xlat(ii)-lft_lat(j)).le.delta)
        yllabels(j) = fabs(lft_lat(j)) + ""
        ylvalues(j) = yndc(ii(0))
        if(lft_lat(j).lt.0) then
          yllabels(j) = yllabels(j) + "~S~o~N~S"
        end if
        if(lft_lat(j).gt.0) then
          yllabels(j) = yllabels(j) + "~S~o~N~N"
        end if
      end if
      delete(ii)
    end do
    bres@tmYLMode   = "Explicit"
    bres@tmYLValues = ylvalues
    bres@tmYLLabels = get_res_value(res2,"tmYLLabels",yllabels)
    bres@tmYLLabelStride = 2
  else
    bres@tmYLOn       = False
    bres@tmYLLabelsOn = False
  end if

;---Right axis tickmarks
  if(isatt(res2,"tmYRValues")) then
    rgt_lat    = get_res_value(res2,"tmYRValues",-1)
    nlat       = dimsizes(rgt_lat)
    yrvalues = new(nlat,float)
    yrlabels = new(nlat,string)

    xfix  = vpx + vpw - 0.0001 ; Just a smidge into the plot to make sure we don't
                               ; get missing values returned.
;
; Loop across each right latitude value that we want a tickmark for,
; and try to find the closest X,Y NDC coordinate pair along this axis.
;
    do j=0,dimsizes(rgt_lat)-1
      NhlNDCToData(plot,xfix,yndc,xlon,xlat)
      ii = minind(fabs(xlat-rgt_lat(j)))
      if(.not.any(ismissing(ii)).and.fabs(xlat(ii)-rgt_lat(j)).le.delta)
        yrlabels(j) = fabs(rgt_lat(j)) + ""
        yrvalues(j) = yndc(ii(0))
        if(rgt_lat(j).lt.0) then
          yrlabels(j) = yrlabels(j) + "~S~o~N~S"
        end if
        if(rgt_lat(j).gt.0) then
          yrlabels(j) = yrlabels(j) + "~S~o~N~N"
        end if
      end if
      delete(ii)
    end do
    bres@tmYROn       = True
    bres@tmYRLabelsOn = True
    bres@tmYUseLeft   = False
    bres@tmYRMode     = "Explicit"
    bres@tmYRValues   = yrvalues
    bres@tmYRLabels   = get_res_value(res2,"tmYRLabels",yrlabels)
  else
    bres@tmYUseLeft   = False
    bres@tmYROn       = False
    bres@tmYRLabelsOn = False
  end if

;---Top axis tickmarks
  if(isatt(res2,"tmXTValues")) then
    top_lon    = get_res_value(res2,"tmXTValues",-1)
    nlon       = dimsizes(top_lon)
    xtvalues = new(nlon,float)
    xtlabels = new(nlon,string)

    yfix  = vpy - 0.0001 ; Just a smidge into the plot to make sure we don't
                         ; get missing values returned.
;
; Loop across each top longitude value that we want a tickmark for,
; and try to find the closest X,Y NDC coordinate pair along this axis.
;
    do j=0,dimsizes(top_lon)-1
      NhlNDCToData(plot,xndc,yfix,xlon,xlat)
      ii = minind(fabs(xlon-top_lon(j)))
      if(.not.any(ismissing(ii)).and.fabs(xlon(ii)-top_lon(j)).le.delta)
        xtlabels(j) = fabs(top_lon(j)) + ""
        xtvalues(j) = xndc(ii(0))
        if(top_lon(j).lt.0) then
          xtlabels(j) = xtlabels(j) + "~S~o~N~W"
        end if
        if(top_lon(j).gt.0) then
          xtlabels(j) = xtlabels(j) + "~S~o~N~E"
        end if
      end if
      delete(ii)
    end do
    bres@tmXTOn       = True
    bres@tmXTLabelsOn = True
    bres@tmXUseBottom = False
    bres@tmXTMode     = "Explicit"
    bres@tmXTValues   = xtvalues
    bres@tmXTLabels   = get_res_value(res2,"tmXTLabels",xtlabels)
  else
    bres@tmXUseBottom = False
    bres@tmXTOn       = False
    bres@tmXTLabelsOn = False
  end if

;---Bottom axis tickmarks
  if(isatt(res2,"tmXBValues")) then
    bot_lon    = get_res_value(res2,"tmXBValues",-1)
    nlon       = dimsizes(bot_lon)
    xbvalues = new(nlon,float)
    xblabels = new(nlon,string)

    yfix  = vpy-vph + 0.0001 ; Just a smidge into the plot to make sure 
                             ; we don't get missing values returned.
;
; Loop across each bottom longitude value that we want a tickmark for,
; and try to find the closest X,Y NDC coordinate pair along this axis.
;
    do j=0,dimsizes(bot_lon)-1
      NhlNDCToData(plot,xndc,yfix,xlon,xlat)
      ii = minind(fabs(xlon-bot_lon(j)))
      if(.not.any(ismissing(ii)).and.fabs(xlon(ii)-bot_lon(j)).le.delta)
        xblabels(j) = fabs(bot_lon(j)) + ""
        xbvalues(j) = xndc(ii(0))
        if(bot_lon(j).lt.0) then
          xblabels(j) = xblabels(j) + "~S~o~N~W"
        end if
        if(bot_lon(j).gt.0) then
          xblabels(j) = xblabels(j) + "~S~o~N~E"
        end if
      end if
      delete(ii)
    end do
    bres@tmXBMode   = "Explicit"
    bres@tmXBValues = xbvalues
    bres@tmXBLabels = get_res_value(res2,"tmXBLabels",xblabels)
    bres@tmXBLabelStride = 2
  else
    bres@tmXBOn       = False
    bres@tmXBLabelsOn = False
  end if

;
; Now that we are done figuring out where to put tickmarks, and
; what labels to use, get any "tm" resources that might have been
; set by the user, and create a blank plot with thes new tickmarks.
;

;---Get rest of user resources that were set with "tm".
  bres = get_res_eq(res2,"tm")

  bres = True   ; Above call will set bres to True if no "tm" resources, so
                ; make sure it is True still.

  bres@gsnDraw  = False
  bres@gsnFrame = False

;
; Create blank plot with new tickmarks (don't use gsn_csm_blank_plot, 
; because it wants to scale the size of your X and Y axes.)
;
  blank  = gsn_blank_plot(wks,bres)

;
; Attach new tickmarks to original plot. This will allow resizing
; if desired. The default is to attach one plot to the center of
; the other one. These two plots are already the same size.
;
  annoid = gsn_add_annotation(plot,blank,False)

;
; Be sure to return the annotation id, otherwise the
; tickmarks will disappear.
;
  anno_str = unique_string("annoid")
  plot@$anno_str$ = annoid

  return(plot)
end
begin

;  files = systemfunc ("ls wrfout_d03_2008*:00")
  files = systemfunc ("ls wrfout_d01_2008-04-08_12:00:00 wrfout_d01_2008-04-08_18:00:00") 

   plot = new(2,graphic)  
   dummy = new(2,graphic)
  figurename = "wrfout_d03_ColumnLiquidIce_Figure2"
  wks = gsn_open_wks("eps" ,figurename)           ; ps,pdf,x11,ncgm,eps
  gsn_define_colormap(wks,"BlAqGrYeOrReVi200") ; select color map

 do ifile = 0, dimsizes(files)-1
  file_target = files(ifile) 
  f     = addfile(file_target+".nc","r")
  fmap  = f

  Height_stag    =(/ (fmap->PH(0,:, :,:)+fmap->PHB(0,:, :,:))/9.81/)
  dims = dimsizes(Height_stag)
  Height_delta= Height_stag(1:dims(0)-1,:,:)-Height_stag(0:dims(0)-2,:,:)
 
  Times_char = f->Times
   Times_char(0,13) = Times_char(0,4) 
   Times_char(0,16) = Times_char(0,4) 
  QCLOUD_column     = f->QCLOUD(0,:,:,:)    ;(Time, bottom_top, south_north, west_east) ;
;  QCLOUD     = (/f->QCLOUD(0,:,:,:)+f->QGRAUP(0,:,:,:)+f->QICE(0,:,:,:)+f->QRAIN(0,:,:,:)+f->QSNOW(0,:,:,:)/)    ;(Time, bottom_top, south_north, west_east) ;
  ALT     = f->ALT(0,:,:,:)
  QCLOUD_column     = (/(f->QCLOUD(0,:,:,:)+f->QGRAUP(0,:,:,:)+f->QICE(0,:,:,:)+f->QSNOW(0,:,:,:))*1000/ALT*Height_delta/) 
  integratedQCLOUD = dim_sum_Wrap( QCLOUD_column(south_north|:, west_east|:, bottom_top|:) )
   v     = f->V10(0,:,:)
   u     = f->U10(0,:,:)

      gsres               = True
      gsres@gsMarkerIndex = 16          ; circle at first
      gsres@gsMarkerThicknessF = 1
      gsres@gsMarkerSizeF = 0.01
      gsres@gsMarkerColor      = "black"


  res                       = True             ; plot mods desired
  res@mpDataBaseVersion  = "Ncarg4_1"
  res@mpDataSetName      = "Earth..4"

  res@tmXBLabelFontHeightF = 0.02
  res@tmYLLabelFontHeightF = 0.02
  res@gsnFrame  = False
  res@gsnDraw  = 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@cnFillMode       = "RasterFill" 
  res@cnLevelSelectionMode  = "ManualLevels"
  res@cnMaxLevelValF        = 200 
  res@cnMinLevelValF       =  0 
  res@cnLevelSpacingF   =  20
  res@lbLabelStride         =2 
  res@lbLabelBarOn        = False           ; turn off individual cb's

  res@tmYLLabelStride = 2
  res@tmXBLabelStride = 2 
  WRF_map_c(fmap, res, 0)                   ; reads info from file
  res@mpOutlineBoundarySets = "AllBoundaries"
   res@vcRefMagnitudeF          = 8.             ; define vector ref mag
   res@vcRefLengthF             = 0.05          ; define length of vec ref
   res@vcRefAnnoOrthogonalPosF  = -1.            ; move ref vector
   res@vcRefAnnoParallelPosF  = 0.999            ; move ref vector
   res@vcMinDistanceF           = 0.04            ; larger means sparser
   res@vcLineArrowHeadMaxSizeF  = 0.0075          ; default: 0.05 (LineArrow), 0.012 (CurlyVector) 
   res@vcGlyphStyle             = "CurlyVector"   ; default: "LineArrow"  
   res@vcRefAnnoFontHeightF   = 0.02 
   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@tmXBOn = False            ; turn off top   labels 
; res@tmYLOn = False            ; turn off right labels
 if (ifile.eq.0) then 
   res@tmXBOn = False            ; turn off right labels
   res@vcRefAnnoOn = True 
  tmres          = True
  tmres@tmYLValues = fspan(70,72,5)
 else
   res@tmXBOn =  True; turn off right labels
   res@vcRefAnnoOn =False 
  tmres@tmXBValues = ispan(-158,-153,1)
 end if 
; res@tiMainString     = "Valid " + chartostring(fmap->Times) + "(UTC) (Init: 2006083018; "+ifile+"hr forecast)"  
 res@gsnLeftStringFontHeightF = 0.025
 res@gsnRightStringFontHeightF = 0.025
; res@gsnLeftString    =   "Condensed water path "+ chartostring(fmap->Times)
 res@gsnLeftString    =  "" ;  chartostring(fmap->Times)
 res@gsnRightString    = "" ;"g m~S~-2" 
  res@cnInfoLabelOn       = False 


  res@cnInfoLabelOrthogonalPosF = -0.04
  res@cnInfoLabelString   = "Min= $ZMN$ Max= $ZMX$"
  lat = f->XLAT(0,:,:)
  longitude = f->XLONG(0,:,:)
   dim_domain=dimsizes(lat)
;  res@mpRightCornerLatF = lat(dim_domain(0)-1,dim_domain(1)-80)
;   res@mpRightCornerLonF = longitude(dim_domain(0)-1,dim_domain(1)-80)

  plot(ifile)                 = gsn_csm_vector_scalar_map(wks,u,v,integratedQCLOUD,res)
;  getvalues plot                         ; retrieve some of the plot resources
;     "tmXBValues"  : tmXBValues          ; values used by NCL at major tick marks
;  end getvalues
  dummy(ifile) = gsn_add_polymarker(wks,plot(ifile),-156.766389,71.295556,gsres)
  tres       =  True
  tres@txFontHeightF = 0.025
  dummy2 = gsn_add_text(wks,plot(ifile),"Barrow",-156.766389,71.295556+0.07,tres)
  tres@txFontHeightF = 0.038
  if(ifile.eq.0) then 
      dummy3 = gsn_add_text(wks,plot(ifile),"a",-158.5,72,tres)
    else 
      dummy3 = gsn_add_text(wks,plot(ifile),"b",-158.5,72,tres)
     ;---Attach the new map tickmarks.
;  print(res)
  end if 
      map = add_map_tickmarks(wks,plot(ifile),tmres)

;  print(tmXBValues)
;  draw(plot(ifile) )
;  frame(wks)
 end do ; ifile

  resP                     = True                ; modify the panel plot
;  resP@txString            = "A plot with a common label bar"
  resP@gsnPanelLabelBar    = True                ; add common colorbar
  resP@lbLabelFontHeightF = 0.02
  resP@lbLabelStride         =2
  resP@lbTitleOn        =  True                ; turn on title
  resP@lbTitleString    = "      g m~S~-2" ;"~C~~V10~g m~S~-2"   ; "O~B~3, ~N~ppbv"             ; title string
  resP@lbTitlePosition  = "Top"              ; title position
  resP@lbTitleFontHeightF= .018                ; make title smaller
  resP@lbTitleDirection = "Across"             ; title direction
  resP@lbLeftMarginF = -0.15
  resP@lbTitleOffsetF = -0.04
  resP@lbLabelOffsetF = 0.02
  resP@lbTitleJust = "CenterCenter"
  resP@lbOrientation = "Vertical"
  resP@gsnPanelBottom = 0.2 
;  resP@lbBoxMajorExtentF = 0.8
   resP@gsnPanelDebug = True
   resP@pmLabelBarWidthF = 0.1
   resP@pmLabelBarHeightF = 0.6
  gsn_panel(wks,plot,(/2,1/),resP)               ; now draw as one plot

 system("convert -trim "+figurename+".eps ../../"+figurename+".png")
 system("mv "+figurename+".eps ../../"+figurename+".eps")
 print("finish plotting "+figurename+".eps")

end