;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; NCL script to plot x-z cross section of radar data
; in netCDF format.
; Use Radx software function RadxConvert to get data 
; in .nc format before proceeding.
; 
; Code copied from Sara Ganetis, 2 May 2016
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

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

xsec_loc = "KDOX"
;casedir = "20150328"
;case = "20150327"

begin
  
  diri = "/media/DATA/rconnelly/Research/"+casedir+"_Case/radar/netcdf_3D/"+xsec_loc+"/"+case+"/"
  fils = systemfunc ("csh -c 'cd " + diri + " ; ls *.nc'")  
  do j = 0,dimsizes(fils)-1,1
       fili = fils(j)    ;+".nc"
       rf    = addfile(diri+fili,"r")
       rd    = rf->REF(:,:,:,:)          ; float REF(time, z0, y0, x0) ; RHO, ZDR
       rd_single_time = rd(0,:,:,:)
       rz    = rf->z0
       rx    = rf->x0
       ry    = rf->y0
       rt    = rf->time
       rlat  = rf->lat0
       rlon  = rf->lon0
       ;print(rlat)
       rlat@_FillValue = max(rlat)
       rlon@_FillValue = rlat@_FillValue
       rd@_FillValue   = 9.96921e+36
       rd@lat2d = rlat          
       rd@lon2d = rlon 
       rz_3d = rd
       ;rz = rz*1000    ; km to m
       rz@units = "km"
       rz_3d!0 = "time"
       rz_3d@units = "km"
       rz_3d@lat2d = rlat
       rz_3d@lon2d = rlon
       rz_3d!0 = "time"

       pl = systemfunc("echo " + fils(j) + " | awk '{print substr($0,5,15)}'")
       year = systemfunc("echo " + fils(j) + " | awk '{print substr($0,5,4)}'")
       month = systemfunc("echo " + fils(j) + " | awk '{print substr($0,9,2)}'")
       day = systemfunc("echo " + fils(j) + " | awk '{print substr($0,11,2)}'")
       hour = systemfunc("echo " + fils(j) + " | awk '{print substr($0,14,2)}'")
       minute = systemfunc("echo " + fils(j) + " | awk '{print substr($0,16,2)}'")
       second = systemfunc("echo " + fils(j) + " | awk '{print substr($0,18,2)}'")
       print("Working on time: " + pl)
       pltName = "xsec_"+xsec_loc+"_"+pl                  ;x_obsdbz, x_obscc, x_obszdr       
       pltType = "png"                             
       pltDir  = "/media/DATA/rconnelly/Research/"+casedir+"_Case/gallery-images/"+xsec_loc+"_xsec/"
       wks = gsn_open_wks(pltType, pltDir+pltName) 

cmap = (/"gray100", "black", "gray100", "lightskyblue1", "lightskyblue", "cornflowerblue", "royalblue3", "blue3", "lawngreen", "green2", "green3", "green4", "darkgreen", "yellow", "gold", "darkgoldenrod1", "darkorange", "darkorange2", "firebrick1", "firebrick3", "firebrick", "firebrick4", "maroon4", "violet"/)


       gsn_define_colormap(wks,cmap)

       ;gsn_define_colormap(wks,"GHRSST_anomaly");    ;for CC/rho-hv

       ; Get height info for labels
       zmin = 0.
       zmax = 10.
       nz   = floattoint(zmax/2 + 1)


if (xsec_loc .eq. "KENX") then

        lonNW = -75.1
        latNW = 43.7

        lonSE = -72.7
        latSE = 41.8

else if (xsec_loc .eq. "KGYX") then

        lonNW = -71.3
        latNW = 44.4

        lonSE = -69.3
        latSE = 42.9

else if (xsec_loc .eq. "KDOX") then

 	lonNW = -75.64
        latNW = 39.4

        lonSE = -75.17
        latSE = 38.24

end if
end if
end if


       ;blat = 39.75
       ;blon = -75.55
       ;bprimlat = 39.2
       ;bprimlon = -72.
       ;alat = 40.65
       ;alon = -75.44
       ;aprimlat = 40.49
       ;aprimlon = -71.86
       lat = (/ latNW, latSE /) ; user specified coordinate pairs
       lon = (/ lonNW, lonSE /)
       nm = getind_latlon2d (rlat,rlon, lat, lon) ; return 2d subscripts
       plane = new(4,float)
       plane = (/ nm(0,1),nm(0,0), nm(1,1),nm(1,0) /) ; approx. start x;y and end x;y point

printVarSummary(rd_single_time)  

       opts = True
       dbz_plane = wrf_user_intrp2d(rd_single_time(:,:,:),plane,0.,opts)
    
;       dbz_plane_1 = wrf_user_intrp3d(rd(0,:,:,:),rz,"v",plane,0.,True)

printVarSummary(dbz_plane)  
;printVarSummary(dbz_plane_1)  

       radres             = True
       radres@gsnMaximize    = True
       radres@cnFillMode    = "RasterFill" 
       radres@gsnSpreadColors = True  
       radres@cnLevelSelectionMode = "ManualLevels"        ; manual contour levels
       radres@trYMinF = 0.
       radres@trYMaxF = 10.
       radres@tmYLMode        = "Explicit"	
       radres@tmYLValues      = (/2.0,4.0,6.0,8.0,10.0/)
       radres@tmYLLabels      = "" + sprintf("%.1f",radres@tmYLValues) 	
       radres@cnMinLevelValF = 10         ;dbz_Dec10 = 0    ; CC = 0.75     ; ZDR = -0.75
       radres@cnMaxLevelValF = 50        ;dbz_Dec10 = 70   ; CC = 1        ; ZDR =  5.5
       radres@cnLevelSpacingF = 2        ;dbz_Dec10 = 5    ; CC = 0.005    ; ZDR =  0.25
       radres@cnLinesOn        = False
       radres@cnLineLabelsOn       = False                 ; turn on line labels
       radres@cnFillOn             = True                  ; turn on color fill
       radres@tiYAxisString      = "Height (km)"
       radres@tiXAxisString     = "Horizontal Distance (km)"
       radres@tiMainString        = xsec_loc+" 2.0 km Reflectivity at " +year+ "-" +month+ "-" +day+ ": " +hour+ ":" +minute+ ":" +second+ " UTC"

       plot =  gsn_csm_contour(wks,dbz_plane(0:98,:),radres)

       ;system("cp xsec_dbz_ilg_"+pl+".png OKX_xDEC10_"+j+".png")
  end do 

end