; This script is used to plot 
; get each hour average value for 24 hours
; get data along two point section
; plot hovmoller for time vs grid along section 

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/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"

begin
; The WRF ARW input file.  
; This needs to have a ".nc" appended, so just do it.
 DATADir = "/geodata/data/private/fgvy024/Build_WRF/WRFV6/run/"
 FILES = systemfunc (" ls -1 " + DATADir + "wrfout_d03_2015-07*") ; contain hourly time stamp
 numFILES = dimsizes(FILES)
 a = addfiles(FILES+".nc","r")
 
 type = "png"                ; map type
 wks = gsn_open_wks(type,"july_hovmoller_hourly_ppt")         ; Create a plot workstation
 rain_exp = wrf_user_getvar(a,"RAINNC",-1)
 rain_con = wrf_user_getvar(a,"RAINC",-1)
 rain = rain_exp + rain_con
 Remp = dimsizes(rain)
 nt = Remp(0)
 nnlat = Remp(1)
 nlon = Remp(2)
; get hourly rainfall
  Ree = rain(1:nt-1,:,:) - rain(0:nt-2,:,:) ; to get hourly precipitation
 ; get new variable with full time step to make it multiplable by 24
  Re = new((/nt, nnlat, nlon/),float)
  Re(1:nt-1,:,:) = Ree(:,:, :) 
  Re(0, :, :) = 0 ; first array will have zero value
  dims = dimsizes(Re)
  pA_lon = 86.2233 ; longitude of point A
  pA_lat = 27.88163333 ; lat of point A
  pB_lon = 86.541725 ; similar for point B
  pB_lat = 27.89950278
    pointA = wrf_user_ll_to_ij(a, pA_lon, pA_lat, True) ; get array point
    pointB = wrf_user_ll_to_ij(a, pB_lon, pB_lat, True)
    xA = pointA(0) - 1
    yA = pointA(1) - 1  
    xB = pointB(0) - 1
    yB = pointB(1) - 1 
  opts = True 
  plane = new(4,float) ; for section along two points
  plane = (/  yA,xA, yB,xB  /)

; interpolate value along the section of A and B
 R_plane = wrf_user_intrp2d(Re,plane,0.,opts)
 Re_dims = dimsizes(R_plane)
 ntim     = Re_dims(0)
 ncross     = Re_dims(1)
 ndays       = ntim/24
 Re3d       = reshape(R_plane,(/24,ndays, ncross/)) ;  hours, days, crosssection
 copy_VarAtts(R_plane,Re3d)
 Re3d_hr_avg = dim_avg_n_Wrap(Re3d,1) ; hourly average (monthly, 0 hours average, 1 hour average, etc)
  res                      = True               ; plot mods desired
  res@cnFillOn             = True               ; turn on color fill
  res@cnFillPalette        = "BlWhRe"           ; set color map
  res@tiMainString         = "Hourly Precipitation"   ; title
  res@cnLevelSelectionMode = "ManualLevels"     ; manual contour levels
  res@cnMinLevelValF       = 0.	        ; min level
  res@cnMaxLevelValF       =  2               ; max level
  res@cnLevelSpacingF      = 0.1               ; contour level spacing
  plot = gsn_csm_hov(wks, Re3d_hr_avg(:,:), res) 
end