;***********************************************
; radar_NWCT.ncl
; Note, need to find slope of kdp line +/- 4 from gate location. Not done yet.
; plots data from *nc files created by the NOAA Weather & Climate Tool Kit
; There are four levels within each variable (different elevation angles).
; The first level is used (index 0, about 0.5 degrees).
; The next index is the azimuth (angle).
; The last index is the range or along radials (gates). 
; the last two maps are based on data that has a simple deattenuation based on
; KDP 
;
; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

begin
  plots_dbz = new ( 6, graphic )
  plots_zdr = new ( 6, graphic )
  plots_cc = new ( 6, graphic )
  plots_dp = new ( 6, graphic )
  plots_dbz_dat = new ( 6, graphic )
  plots_zdr_dat = new ( 6, graphic )
; dir = "./"
  do i_file = 0,5
  filenames = ((/"KOKX_V06_20130209_000138.nc","KOKX_V06_20130209_003053.nc", \
                "KOKX_V06_20130209_010213.nc", "KOKX_V06_20130209_013125.nc", \
                "KOKX_V06_20130209_020046.nc", "KOKX_V06_20130209_023001.nc"/))
  
; map Reflectivity  
  filename = filenames(i_file)
   if (i_file.eq.0)then
   wks = gsn_open_wks("png",filename)           ; send graphics to PNG file
   end if
  print("filename = " + filename)
; varname  = "Reflectivity_HI"
  f        = addfile(filename,"r")
; end do
  var     := f->Reflectivity_HI(0,:,:)
  dimsvar = dimsizes(var)
  print("dimsvar = " + dimsvar)
  scale_factor = var@scale_factor_metadata
  add_offset = var@add_offset_metadata
  dbz      := byte2flt( (/ var /) ) ;; strip data and don’t apply offset / scale yet. 
  copy_VarAtts(var,dbz)
  dbz = where(dbz.lt.0, dbz+256,dbz)
  dbz = (dbz * scale_factor)+add_offset	;; now scale the data
  angles          = f->azimuthR_HI(0,:)
  dims_ang = dimsizes(angles)
  j_s = 0
  j_e = dims_ang(0)-1

  dsizes_dbz       = dimsizes(dbz)
  i_s = 0
  i_e = 4*dsizes_dbz(0)/5
  dbz_sub = dbz(j_s:j_e,i_s:i_e)
  dbz_o = dbz
;---To see a list of all the variables on the file
   if (i_file .eq. 0)then
   print(getfilevarnames(f))
   end if
;
 angles=90-angles
 angles_sub = angles(j_s:j_e)
  DEGTORAD = 0.017453292519943
  xcenter  = 0.0
  ycenter  = 0.0
; Create 2D coordinate arrays.
  x=f->distanceR_HI
  x_sub = x(i_s:i_e)
  angles2d = conform(dbz_sub,angles_sub,0)
  x2d      = conform(dbz_sub,x_sub,1)/1000 ; convert to km
  xarr     =  xcenter + x2d  * cos(DEGTORAD * angles2d)
  yarr     =  ycenter + x2d  * sin(DEGTORAD * angles2d)

; wks = gsn_open_wks("png",str_join((/filename,varname/),"_"))           ; send graphics to PNG file

  dbres                   = True
  dbres@gsnDraw = False
  dbres@gsnFrame = False
  dbres@cnInfoLabelOn = False 

  dbres@lbLabelBarOn = False 
  dbres@gsnMaximize       = True
  dbres@sfXArray          = xarr
  dbres@sfYArray          = yarr
  dbres@cnFillOn          = True
  dbres@cnFillPalette     = "NCV_jet"         ; set color map
  dbres@cnLinesOn         = False
  dbres@cnFillMode        = "RasterFill"      ; this mode is fastest
  dbres@trGridType        = "TriangularMesh"
  dbres@tiMainString      = filename
  dbres@tiXAxisString     = "Zonal Distance (km)"
  dbres@tiYAxisString     = "Meridional Distance (km)"
  dbres@lbOrientation     = "Vertical"
  dbres@cnLevelSelectionMode = "ExplicitLevels" ; use explicit levels
  dbres@cnLevels             = ispan(-10,55,5) ; set the contour levels
  dbres@cnLineLabelsOn       = False            ; do not use line labels
  dbres@cnFillOn             = True             ; color fill
  dbres@cnLinesOn            = False            ; do not draw contour lines
  contour = gsn_csm_contour(wks,dbz_sub,dbres)
  plots_dbz(i_file) = contour
; will reuse old dbz variable for differential reflectivity
  delete(dbz)
  delete(var)
  var      = f->DifferentialReflectivity_HI(0,:,:)
  scale_factor = var@scale_factor_metadata
  add_offset = var@add_offset_metadata
  dbz      = byte2flt( (/ var /) ) ;; strip data and don’t apply offset / scale yet.
  copy_VarAtts(var,dbz)
  dbz = where(dbz.lt.0, dbz+256,dbz)
  dbz = (dbz * scale_factor)+add_offset	;; now scale the data
  dbz_sub = dbz(j_s:j_e,i_s:i_e)
  zdr_o = dbz
  angles2d = conform(dbz_sub,angles_sub,0)
  x2d      = conform(dbz_sub,x_sub,1)/1000 ; convert to km
  dfres                   = True
  dfres@lbLabelBarOn = False 
  dfres@gsnDraw = False
  dfres@gsnFrame = False
  dfres@cnInfoLabelOn = False 
  dfres@gsnMaximize       = True
  dfres@sfXArray          = xarr
  dfres@sfYArray          = yarr
  dfres@cnFillOn          = True
  dfres@cnFillPalette     = "NCV_jet"         ; set color map
  dfres@cnLinesOn         = False
  dfres@cnFillMode        = "RasterFill"      ; this mode is fastest
  dfres@trGridType        = "TriangularMesh"
  dfres@tiMainString      = filename
  dfres@tiXAxisString     = "Zonal Distance (km)"
  dfres@tiYAxisString     = "Meridional Distance (km)"
  dfres@lbOrientation     = "Vertical"
; cnres@tiMainString      = str_join((/filename," - ",varname/),"_")
  dfres@tiMainString      = filename
  dfres@cnLevelSelectionMode = "ExplicitLevels" ; use explicit levels
; dfres@cnLevels             = ispan(-2,6,1) ; set the contour levels
  dfres@cnLevels             = ispan(-2,5,1) ; set the contour levels
  dfres@cnLineLabelsOn       = False            ; do not use line labels
  dfres@cnFillOn             = True             ; color fill
  dfres@cnLinesOn            = False            ; do not draw contour lines
; cnres@gsnSpreadColorStart= -1

  contour = gsn_csm_contour(wks,dbz_sub,dfres)
  plots_zdr(i_file) = contour
; will reuse old dbz variable for correlation coefficient
  delete(dbz)
  delete(var)
  var      = f->CorrelationCoefficient_HI(0,:,:)
  scale_factor = var@scale_factor_metadata
  add_offset = var@add_offset_metadata
;  dbz      = byte2flt( (/ f->Reflectivity_HI(0,:,:) /) ) ;; strip data and don’t apply offset / scale yet.
  dbz      = byte2flt( (/ var /) ) ;; strip data and don’t apply offset / scale yet.
  copy_VarAtts(var,dbz)
  dbz = where(dbz.lt.0, dbz+256,dbz)
  dbz = (dbz * scale_factor)+add_offset	;; now scale the data
  dbz_sub = dbz(j_s:j_e,i_s:i_e)
  angles2d = conform(dbz_sub,angles_sub,0)
  x2d      = conform(dbz_sub,x_sub,1)/1000 ; convert to km
  ccres                   = True
  ccres@lbLabelBarOn = False 
  ccres@gsnDraw = False
  ccres@gsnFrame = False
  ccres@cnInfoLabelOn = False 
  ccres@gsnMaximize       = True
  ccres@sfXArray          = xarr
  ccres@sfYArray          = yarr
  ccres@cnFillOn          = True
  ccres@cnFillPalette     = "NCV_jet"         ; set color map
  ccres@cnLinesOn         = False
  ccres@cnFillMode        = "RasterFill"      ; this mode is fastest
  ccres@trGridType        = "TriangularMesh"
  ccres@tiMainString      = filename
  ccres@tiXAxisString     = "Zonal Distance (km)"
  ccres@tiYAxisString     = "Meridional Distance (km)"
  ccres@lbOrientation     = "Vertical"
  ccres@tiMainString      = filename
  ccres@cnLevelSelectionMode = "ExplicitLevels"    ; set explicit contour levels
  ccres@cnLevels    = (/ 0.70,0.71,0.72,0.73,0.74,0.75,0.76,0.77,0.78,0.79,0.80,0.81,\
                          0.82,0.83,0.84,0.85,0.86,0.87,0.88,\
                          0.89,0.90,0.91,0.92,0.93,0.94,0.95,0.96,\
                          0.97,0.98,0.99,1.0/)

  ccres@cnFillColors = (/ 21,25,31,35,39,43,48,53,\
                           61,81,85,89,95,99,104,108,115,123,\
                           130,135,143,148,155,159,167,171,\
                           180,184,188,192,196,200,205/)

  ccres@cnLineLabelsOn       = False            ; do not use line labels
  ccres@cnFillOn             = True             ; color fill
  ccres@cnLinesOn            = False            ; do not draw contour lines
; cnres@gsnSpreadColorStart= -1

; x = ispan (1,10,1)
  contour = gsn_csm_contour(wks,dbz_sub,ccres)
  plots_cc(i_file) = contour
; will reuse old dbz variable for Differential Phase
  delete(dbz)
  delete(var)
  var      = f->DifferentialPhase_HI(0,:,:)
  scale_factor = var@scale_factor_metadata
  add_offset = var@add_offset_metadata
;  dbz      = byte2flt( (/ f->Reflectivity_HI(0,:,:) /) ) ;; strip data and don’t apply offset / scale yet.
  dbz      = byte2flt( (/ var /) ) ;; strip data and don’t apply offset / scale yet.
  copy_VarAtts(var,dbz)
  dbz = where(dbz.lt.0, dbz+256,dbz)
  dbz = (dbz * scale_factor)+add_offset	;; now scale the data
  wk_smooth121(dbz)
  dbz_sub = dbz(j_s:j_e,i_s:i_e)
  dfp_o = dbz
  angles2d = conform(dbz_sub,angles_sub,0)
  x2d      = conform(dbz_sub,x_sub,1)/1000 ; convert to km
  dpres                   = True
  dpres@lbLabelBarOn = False 
  dpres@gsnDraw = False
  dpres@gsnFrame = False
  dpres@cnInfoLabelOn = False 
  dpres@gsnMaximize       = True
  dpres@sfXArray          = xarr
  dpres@sfYArray          = yarr
  dpres@cnFillOn          = True
  dpres@cnFillPalette     = "NCV_jet"         ; set color map
  dpres@cnLinesOn         = False
  dpres@cnFillMode        = "RasterFill"      ; this mode is fastest
  dpres@trGridType        = "TriangularMesh"
  dpres@tiMainString      = filename
  dpres@tiXAxisString     = "Zonal Distance (km)"
  dpres@tiYAxisString     = "Meridional Distance (km)"
  dpres@lbOrientation     = "Vertical"
  dpres@tiMainString      = filename
  dpres@cnLevelSelectionMode = "ExplicitLevels"    ; set explicit contour levels
; dfres@cnLevels             = ispan(15,180,15) ; set the contour levels
; ccres@cnLevels    = (/ 0.75,0.76,0.77,0.78,0.79,0.80,0.81,\
;                         0.82,0.83,0.84,0.85,0.86,0.87,0.88,\
;                         0.89,0.90,0.91,0.92,0.93,0.94,0.95,0.96,\
;                         0.97,0.98,0.99,1.0,1.01,1.02,1.03,1.04,1.05/)

  dpres@cnLevels    = (/ 10,20,30,40,50,60,70,80,90,120,180/)
; ccres@cnFillColors = (/ 21,25,31,35,39,43,48,53,\
;                          61,81,85,89,95,99,104,108,115,123,\
;                          130,135,143,148,155,159,167,171,\
;                          180,184,188,192,196,200,205/)

  dpres@cnLineLabelsOn       = False            ; do not use line labels
  dpres@cnFillOn             = True             ; color fill
  dpres@cnLinesOn            = False            ; do not draw contour lines
; cnres@gsnSpreadColorStart= -1

; x = ispan (1,10,1)
  contour = gsn_csm_contour(wks,dbz_sub,dpres)
  plots_dp(i_file) = contour


; var --> calculate deattenuation
;Bringi et al. (1990): 
;Bringi, V.N., V. Chandrasekar, N. Balakrishnan, and D.S. Zrnić, 1990: An Examination of Propagation Effects in Rainfall on Radar Measurements at Microwave Frequencies. J. Atmos. Oceanic Technol., 7, 829–840, https://doi.org/10.1175/1520-0426(1990)007<0829:AEOPEI>2.0.CO;2
 
  dims_o = dimsizes(dfp_o)
; Calculate KDP (Specific Differential Phase) from Differential Phase.
; the grid spacing is 0.1 km
  dfp_dif = dfp_o
  dfp_dif = 0
; need to relabel kdp with correct units, variable name
; not done here (KDP is specific differential attenuation and has
; units of db/km).

; move away from ground cluster (2 km) 
  do j = 0,dims_o(0)-1
; do j = 0,0,1
  do i = 1,dims_o(1)-1
   dfp_dif(j,i) = (dfp_o(j,i) - dfp_o(0,20))
  end do
  end do
; from Snyder et al (2010): http://journals.ametsoc.org/doi/full/10.1175/2010JTECHA1356.1
; alpha  = 0.313 and alpha_dp 0.0483
  alpha = 0.16
  alpha_dp = 0.36

; no integration along radial in this simple approach
; sum_att = dbz_o
; sum_att = 0.
; do j = 0,dims_o(0)-1
; do i = 1,dims_o(1)-1
;  sum_att(j,i) = sum_att(j,i-1) + alpha*kdp(j,i)
; end do
; end do
; print("calculated PIA")
  dbz_new = dbz_o
  zdr_new = zdr_o
  added = 0
  do j = 0,dims_o(0)-1
  do i = 1,dims_o(1)-1
;  dbz_new(j,i) = dbz_o(j,i) + sum_att(j,i)
   check = ismissing(dbz_o(j,i))
;  if (check.eq.False)then
   if (dfp_dif(j,i)*alpha.gt.0)then
   dbz_new(j,i) = dbz_o(j,i) + alpha*dfp_dif(j,i)
   zdr_new(j,i) = zdr_o(j,i) + alpha_dp*alpha*dfp_dif(j,i)
   added = added + 1
   else
   dbz_new(j,i) = dbz_o(j,i)
   zdr_new(j,i) = zdr_o(j,i)
   end if
;  end if
  end do
  end do
  print("added grid numbers = " + added)
  dbz_sub = dbz_new(j_s:j_e,i_s:i_e)
  contour = gsn_csm_contour(wks,dbz_sub,dbres)
  plots_dbz_dat(i_file) = contour
  zdr_sub = zdr_new(j_s:j_e,i_s:i_e)
  contour = gsn_csm_contour(wks,zdr_sub,dfres)
  plots_zdr_dat(i_file) = contour
  ; Panel the WRF plots.
    pnlres                            = True
;   pnlres@txString                   = t@description + "  (" + t@units + ")"
    pnlres@txString                   = "Observed Reflectivity (dbZ)"
    pnlres@gsnPanelYWhiteSpacePercent = 13       ; Add white space b/w plots.
    pnlres@gsnPanelLabelBar           = True    ; Turn on common labelbar
    pnlres@lbLabelAutoStride          = True    ; Spacing of lbar labels.
    pnlres@lbBoxMinorExtentF          = 0.13
    pnlres@lbLabelFontHeightF    = .015
    pnlres@lbPerimOn             = False
    pnlres@lbLabelFont      = "Helvetica"         ; label font
    pnlres@lbTitleOn        =  True               ; turn on title
    pnlres@lbTitleString    = "dbZ"        ; title string
;   pnlres@lbTitlePosition  = "Top"               ; title position
    pnlres@lbTitlePosition  = "Bottom"               ; title position
    pnlres@lbTitleFontHeightF= .015               ; make title smaller
    pnlres@lbTitleDirection = "Across"            ; title direction
    pnlres@pmLabelBarWidthF      = 0.5
    pnlres@pmLabelBarHeightF     = 0.15
    pnlres@lbTopMarginF = 0.001
    pnlres@pmLabelBarOrthogonalPosF = 0.02   ; position wrt plot

    pnlres@lbTitleOffsetF =   -0.4
   end do
    gsn_panel(wks,(/plots_dbz/),(/2,3/),pnlres)
    pnlres@txString                   = "Observed ZDR (dbZ)"
    gsn_panel(wks,(/plots_zdr/),(/2,3/),pnlres)
    pnlres@txString                   = "Observed Correlation Coeficient"
    gsn_panel(wks,(/plots_cc/),(/2,3/),pnlres)
    pnlres@txString                   = "Observed Differential Phase (Deg)"
    gsn_panel(wks,(/plots_dp/),(/2,3/),pnlres)
    pnlres@txString                   = "Observed Corrected Reflectivity (dbZ)"
    gsn_panel(wks,(/plots_dbz_dat/),(/2,3/),pnlres)
    pnlres@txString                   = "Observed Corrected ZDR (dbZ)"
    gsn_panel(wks,(/plots_zdr_dat/),(/2,3/),pnlres)

  end