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
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
;                    printVarSummary(RF_short)
;                    printMinMax(RF_short,0)      
;                    asciiwrite("RF_short.txt",RF_short)

   RF_flt            =  tofloat(RF_short)
   RF_flt@_FillValue = -9999
                    printVarSummary(RF_flt)
                    printMinMax(RF_flt,True)
                    asciiwrite("RF_flt.txt",RF_flt)
;
    RF_with_missing            = RF_flt/(RF_short@scale_factor)
    RF_with_missing@_FillValue = -9999
                    printVarSummary(RF_with_missing)
                    printMinMax(RF_with_missing,True)
                    asciiwrite("RF_with_missing.txt",RF_with_missing)

    RF_fnl          =  where(RF_with_missing.lt.(-88),RF_flt@_FillValue,RF_with_missing)
;                       copy_VarMeta(RF_flt,RF_fnl)
                     printVarSummary(RF_fnl)
;                     printMinMax(RF_fnl,0)
        do ii      = 0,79
                    print(ii)
                    printMinMax(RF_fnl(:,:,ii),True)
        end do
   RF_fnl@lon2d      =  longitude
   RF_fnl@lat2d      =  latitude
                    asciiwrite("RF_fnl.txt",RF_fnl)

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

  res         =   True ; plot mods desired
  res@cnFillOn =  True ; enable contour fill
  res@cnConstFEnableFill =  True ; enable contour fill
  res@cnLevelSelectionMode =  "ExplicitLevels"     
  res@cnLevels             =  (/2,3,4,5,6,10,12/)
  res@mpMinLatF            = 0.0
  res@mpMaxLatF            = 25.0
  res@mpMinLonF            = 60.0
  res@mpMaxLonF            = 100.0
  res@tiMainString         = file_name
  res@gsnRightString       = "Layer 61"
;    plot  =    gsn_contour(XWKS,RF_fnl(:,:,0),res)
    plot  =    gsn_csm_contour_map(XWKS,RF_fnl(:,:,61),res)
end