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

    initdate = "2011_02_15_06"

    wrfdir = "/psd2QPF/egrelldata/HMT/wrfout/noCuPa_v3.6/2011_02_15_06/arw-mor/"
    filnam = "wrfout_d01_2011-02-15_18:30:00"
    pltdir = wrfdir

    zmin = 0.
    zmax = 8.           ; Max extent of vertical axis, in km
    nz   = floattoint(zmax + 1)

    ixspts = (/134,163,202,194/)    ; pts for cross section (xA,yA,xB,yB)
    xspts = 1.0*ixspts
;=================================================================
;=================================================================
; read wrfout file

    print("Reading file: "+filnam)
    wrf = addfile(wrfdir+filnam+".nc","r")     ; load the requested file

    lon = wrf->XLONG(0,:,:)                  ; longitude

; first read in z, and determine the vertical levels needed to 
; get to max requested height.  Then only read in the levels needed.
; (This is an attempt to get a smoother lower boundary by getting
; higher vertical resolution, since the # levels in the interpolated 
; field is fixed.)

    z1 = wrf_user_getvar(wrf,"z",0)             ; 3d Z
    zmaxm = zmax*1000.
    
    dimz = dimsizes(z1)
    kmax = dimz(0)       ; grid dimension in vertical direction
    ksave = 1
    do k=2,kmax-1
      if (any(z1(k,:,:).lt.zmaxm)) then
        ksave = k
      end if
    end do 
    ksave = ksave+2    ; ensure that we have enough data for interpolation
    print("kmax = "+kmax+"      ksave = "+ksave)

    z = z1(0:ksave,:,:)

    print("max(z(0,:,:)) = "+max(z(0,:,:)))
    print("min(z(0,:,:)) = "+min(z(0,:,:)))

    printVarSummary(z)
    delete(z1)


    refl = wrf->REFL_10CM(0,0:ksave,:,:)          ; model reflectivity

    print("Completed read of WRF file")

;========================================================================
; interpolate to plane of cross section

    print("xspts = "+xspts)
    dbz_plane = wrf_user_intrp3d(refl,z,"v",xspts,0.,True)
    zz        = wrf_user_intrp3d(z,z,"v",xspts,0.,True)

    print("max(zz(0,:)) = "+max(zz(0,:)))
    print("min(zz(0,:)) = "+min(zz(0,:)))
    print("zz@_FillValue = "+zz@_FillValue)
    print("max(zz) = "+max(zz))
    print("min(zz) = "+min(zz))

    printVarSummary(zz)
    printVarSummary(dbz_plane)

    print("max(dbz_plane(0,:)) = "+max(dbz_plane(0,:)))
    print("min(dbz_plane(0,:)) = "+min(dbz_plane(0,:)))
    print("dbz_plane@_FillValue = "+dbz_plane@_FillValue)
    print("max(dbz_plane) = "+max(dbz_plane))
    print("min(dbz_plane) = "+min(dbz_plane))

    X_plane = wrf_user_intrp2d(lon,xspts,0.,True)

    bb = ind(zz(:,0) .gt. zmax*1000.)    ; find index of first pt over 8000m
    zlim = bb(0) - 1
    print("zlim = "+zlim)

; X-axis 
    dimsX = dimsizes(X_plane)
    xspan = dimsX(0)-1
    print("xspan = "+xspan)

;========================================================================
; Begin plot info section
;========================================================================
; set generic resources for all plots

     res = True

     res@tiXAxisString           = " "
     res@tiYAxisString           = "Height (km)"

     res@cnMissingValPerimOn     = True
     res@cnMissingValFillColor   = "gray"
     res@cnMissingValFillPattern = 0
     res@cnLineLabelBackgroundColor = -1

     res@tmXTOn                  = False
     res@tmYROn                  = False
     res@tmXBMode                = "Explicit"
     res@tmXBValues              = fspan(0,xspan,2)           ; Create tick marks
     posA = "(" + sprinti("%3.0i",ixspts(0)) + "," + sprinti("%3.0i",ixspts(1)) + ")"
     posB = "(" + sprinti("%3.0i",ixspts(2))+","+sprinti("%3.0i",ixspts(3)) + ")"
     res@tmXBLabels              = (/posA,posB/)
     res@tmYLMode                = "Explicit"
     res@tmYLValues              = fspan(0,zlim,nz)             ; Create tick marks
     res@tmYLLabels              = sprintf("%.1f",fspan(zmin,zmax,nz)) 

     res@gsnRightString = "Stn: BBY"


; resources for DBZ
   resDBZ = res
   resDBZ@gsnLeftString = "REFL_10cm (dBZ)"
   resDBZ@cnFillOn        = True
   resDBZ@cnLinesOn       = False
   resDBZ@cnLineLabelsOn  = False
   resDBZ@cnInfoLabelOn   = False     ; No info box "contour from... to.."
   resDBZ@gsnSpreadColors = True
   resDBZ@cnLevelSelectionMode = "ManualLevels"
   resDBZ@cnMinLevelValF  = -35.
   resDBZ@cnMaxLevelValF  =  40.
   resDBZ@cnLevelSpacingF =  1.0
;!!!!!!!!!!  commenting out the following line will give a 
; different fill value for missing/underground data
   resDBZ@cnFillMode      = "RasterFill"

;========================================================================
;========================================================================

; plot reflectivity

    wks = gsn_open_wks("png",pltdir+"/xs_refl_test_raster")
    gsn_define_colormap(wks,"BlAqGrYeOrRe") 

    contour_refl = gsn_csm_contour(wks,dbz_plane(0:zlim,:),resDBZ)

 end