load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"

  aDataDir       = "./"
  aFileName      = "wrfout_d01_2003-07-15_00:00:00.nc"
  aFile          = aDataDir+aFileName
  a = addfile(aFile,"r")                            ; a is the wrf output with the format of rcm

  bDataDir       = "/Users/shea/Data/TRMM/3B42RT/"
  bFileName      = "3B42RT.2007062518.nc"
  bFile          = bDataDir+bFileName
  b = addfile(bFile,"r")                            ; b is the trmm data with the format of grid
; --- Read WRF coordinates
  lat2d = wrf_user_getvar(a,"XLAT",0)               ; 	[south_north] x [west_east]
  lon2d = wrf_user_getvar(a,"XLONG",0)
  printVarSummary(lat2d)

; --- TRMM lon are 0-360; make WRF 0-360
  lon2d = where(lon2d.lt.0, lon2d+360, lon2d) 


; --- WRF boundaries
  minLatWrf  = min(lat2d)      
  maxLatWrf  = max(lat2d)
  minLonWrf  = min(lon2d)
  maxLonWrf  = max(lon2d)

; --- get 'extra' grid points; extra is user specified [ >=0 ]
  extra      = 2            
  minLatWrfx = minLatWrf - extra
  maxLatWrfx = maxLatWrf + extra
  minLonWrfx = minLonWrf - extra
  maxLonWrfx = maxLonWrf + extra

; --- Read TRMM data for region  
  prcTrmm   = short2flt(b->PRC(:,{minLatWrfx:maxLatWrfx},{minLonWrfx:maxLonWrfx}))
  printVarSummary(prcTrmm)   ; <====== **LOOK** at this!
  print("prcTrmm: min="+min(prcTrmm)+"   max="+max(prcTrmm) )

; --- Regrid
  prcWrf    = rgrid2rcm_Wrap(prcTrmm&lat, prcTrmm&lon, prcTrmm, lat2d, lon2d,0)
  prcWrf@long_name = "TRMM on WRF grid"
  printVarSummary(prcWrf)
  print("prcWrf: min="+min(prcWrf)+"   max="+max(prcWrf) )

; +++ PLOT
  ;dimWrf  = dimsizes(lat2d)
  ;nlat    = dimWrf(0)
  ;mlon    = dimWrf(1)

   wksType = "png"
   wksName ="NCL_rain_2rcm"
   wks     = gsn_open_wks(wksType,wksName)

   res                     = True
   res@gsnMaximize         = True
   res@gsnAddCyclic        = False
   prcWrf@lat2d            = lat2d
   prcWrf@lon2d            = lon2d

   res@mpMinLatF            = minLatWrf
   res@mpMaxLatF            = maxLatWrf
   res@mpMinLonF            = minLonWrf
   res@mpMaxLonF            = maxLonWrf
   res@mpCenterLonF         = 0.5*(res@mpMinLonF+res@mpMaxLonF)
   res@mpFillOn             = False

   res@cnFillOn             = True     ; turn on color
   res@cnFillMode           = "RasterFill"
   res@cnLinesOn            = False    ; turn off contour lines
   res@cnLineLabelsOn       = False    ; turn off contour line labels

   res@cnLevelSelectionMode = "ExplicitLevels"              
   res@cnLevels             = (/0.1,1,2.5,5,10,15,20,25,50,75/) ;   ; 10 contour values
   res@cnFillPalette        = (/"Snow","PaleTurquoise","PaleGreen"\ ; 11 contour colors
                              ,"SeaGreen3" ,"Yellow","Orange"    \
                              ,"HotPink","Orange","HotPink","Red"\
                              ,"Violet", "Purple", "Brown"       /)       

   plot = gsn_csm_contour_map(wks,prcWrf(0,:,:),res)

 ;;print(prcWrf)