;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; 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"

site="KGYX"

begin

  diri = "/media/DATA/rconnelly/Research/"+casedir+"_Case/radar/netcdf_3D/"+site+"/"+case+"/"
  fils = systemfunc ("csh -c 'cd " + diri + " ; ls *.nc'")  
  do j = 0,dimsizes(fils)-1,1
  ; j = 3
       fili = fils(j)    ;+".nc"
       rf    = addfile(diri+fili,"r")
       rd    = rf->REF(:,:,:,:)          ; float REF(time, z0, y0, x0) ; RHO, ZDR
       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 = pl                  ;x_obsdbz, x_obscc, x_obszdr       
       pltType = "png"                             
       pltDir  = "/home/rconnelly/Research/"+casedir+"_Case/gallery-images/"+site+"_2km/"
       wks = gsn_open_wks(pltType, pltDir+"rainbow_"+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,"NCV_bright");    ;for CC/rho-hv
       
; Map plotting options
       pltres = True
       radres = True
       radres@mpFillOn = False
       radres@mpGeophysicalLineColor = "Black"
       radres@mpNationalLineColor    = "Black"
       radres@mpUSStateLineColor     = "Black"
       radres@mpGridLineColor        = 0
       radres@mpLimbLineColor        = "Black"
       radres@mpPerimLineColor       = "Black"
       radres@mpUSStateLineThicknessF = 4.0
       radres@mpGeophysicalLineThicknessF = 4.0
       radres@mpNationalLineThicknessF = 6.0
       radres@mpOutlineBoundarySets = "GeophysicalAndUSStates"
       radres@mpDataSetName = "Earth..4"
       radres@sfXArray = rlon
       radres@sfYArray = rlat
       radres@mpMinLonF = -72.4 ;min(rlon)
       radres@mpMaxLonF = -68. ;max(rlon)-1.8
       radres@mpMinLatF = 42.2; ;min(rlat)+0.5
       radres@mpMaxLatF = 45.5 ;(rlat)-1.3
       radres@mpDataBaseVersion    = "MediumRes"
       radres@gsnMajorLatSpacing = 1
       radres@gsnMajorLonSpacing = 1


; Radar plotting options
       radres             = True
       radres@gsnAddCyclic = False
       radres@gsnMaximize    = True
       radres@cnFillMode    = "RasterFill" 
       radres@gsnSpreadColors = True  
       radres@cnLevelSelectionMode = "ManualLevels"        ; manual contour levels
       radres@cnMinLevelValF = 0         ;dbz_Dec10 = 0    ; CC = 0.75     ; ZDR = -0.75
       radres@cnMaxLevelValF = 60        ;dbz_Dec10 = 70   ; CC = 1        ; ZDR =  5.5
       radres@cnLevelSpacingF = 0.5        ;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@tiMainString        = site+" 2.0 km Reflectivity at " +year+ "-" +month+ "-" +day+ ": " +hour+ ":" +minute+ ":" +second+ " UTC"
       radres@tiMainFontHeightF = 0.02	; default is 0.025
       radres@tiMainFontThicknessF = 0.5	; default is 1.0

       rd_plane = rd(0,2,:,:)	; z is every 0.5 km starting at 1.0 km, so index 2 is 2 km

       rd_plane@lat2d = rlat
       rd_plane@lon2d = rlon

       plot =  gsn_csm_contour_map(wks,rd_plane,radres)

; Construct a polyline to show cross-section location

  gres                 = True
  gres@gsFillColor     = "black"
  gres@cnLineDrawOrder      = "PreDraw"
  gres@gsLineThicknessF     = 5.0

  lat = (/44.4,  42.9/)
  lon = (/-71.3, -69.3/)
  dum = gsn_add_polyline(wks,plot,lon,lat,gres)

  draw(plot)

  end do 

end