load "$NCARG_ROOT/lib/ncarg/nclex/gsun/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
begin

; Read file. 
file_name = "2A25.20090515.65505.7.HDF"
hdf4_file=addfile(file_name, "r") 

; Print metadata information.
;print(hdf4_file)

; Read lat/lon data.
longitude=hdf4_file->Longitude                   ; [nscan | 9247] x [nray | 49] 
latitude=hdf4_file->Latitude
         printMinMax(longitude,True)
         printMinMax(latitude,True)
         print(dimsizes(longitude))
         print(dimsizes(latitude))

; Correct units to follow the CF conventions.
; In the HDF4 file, the attribute values are "degrees". 

longitude@units = "degrees_east"
latitude@units  = "degrees_north"

; Read data for plot.
   RF_short =  hdf4_file->rain                   ; [nscan | 9247] x [nray | 49] x [ncell1 | 80]
                     printVarSummary(RF_short)
                     printMinMax(RF_short,0)      
;                    asciiwrite("RF_short.txt",RF_short)

   RF_dim   = dimsizes(RF_short)
   print(RF_dim)
   nscan    = RF_dim(0)                            ; 9247
   nray     = RF_dim(1)                            ; 49
   ncell    = RF_dim(2)                            ; 80

   RF_flt   = RF_short/tofloat(RF_short@scale_factor)
   RF_flt@_FillValue = -999.0                      
   RF_flt@units      = RF_short@units
   RF_flt@long_name  = RF_short@hdf_name
   copy_VarCoords(RF_short, RF_flt)
                     printVarSummary(RF_flt)
                     printMinMax(RF_flt, True)  

   RF_flt   = where (RF_flt.lt.0, RF_flt@_FillValue, RF_flt)
                    printVarSummary(RF_flt)
                    printMinMax(RF_flt,True)
;                   asciiwrite("RF_with_missing.txt",RF_with_missing)
 

   do ii= 0,ncell-1
      print("ii="+sprinti("%0.2i", ii)+"  RF_flt: min="+min(RF_flt(:,:,ii)) \
                                               +" max="+max(RF_flt(:,:,ii)) )
   end do

;                   asciiwrite("RF_flt.txt",RF_flt)

  XWKS   =   gsn_open_wks("png","Vert_Plot")
             gsn_define_colormap(XWKS,"default")

  RF_flt@lon2d      =  longitude
  RF_flt@lat2d      =  latitude

  res               =   True ; plot mods desired
  res@cnFillOn      =  True ; enable contour fill
;;res@cnFillMode    = "CellFill"    
  res@cnLinesOn     =  False             ; Turn off contour lines
  res@cnConstFEnableFill = True          ; enable contour fill
;  res@cnLevelSelectionMode =  "ExplicitLevels"     
;  res@cnLevels             =  (/2,3,4,5,6,10,12/) ;  res@mpMinLatF            = 0.0
  res@tiMainString         = file_name

    print("----------------")
    nscan_index = 12
    print("Vertical profile: nscan_index="+nscan_index)
    res@gsnRightString = "nscan="+nscan_index
    printMinMax(RF_flt(nscan_index,:,:),0)
    plot  = gsn_contour(XWKS,RF_flt(nscan_index,:,:),res)

    print("----------------")
    print("Swath")
    cell_index = 61
    print("Swatth: cell_index="+cell_index)
    res@gsnRightString = "Layer="+cell_index
    res@mpMinLatF            =  0.0
    res@mpMaxLatF            = 25.0
    res@mpMinLonF            = 60.0
    res@mpMaxLonF            = 100.0
    printMinMax(RF_flt(:,:,cell_index),0)
    plot  = gsn_csm_contour_map(XWKS,RF_flt(:,:,cell_index),res)
end