; Copyright (C) 2014 by The HDF Group.
;   All rights reserved.
;
;  This example code illustrates how to access and visualize GESDISC TRMM
; version 7 HDF4 Level 2 2A12 file in NCL. 
;
;  If you have any questions, suggestions, comments on this example, please 
; use the HDF-EOS Forum (http://hdfeos.org/forums). 
;
;  If you would like to see an  example of any other NASA HDF/HDF-EOS data 
; product that is not listed in the HDF-EOS Comprehensive Examples page 
; (http://hdfeos.org/zoo), feel free to contact us at eoshelp@hdfgroup.org or 
; post it at the HDF-EOS Forum (http://hdfeos.org/forums).

; Tested under: NCL 6.1.2
; Last updated: 2014-4-1

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 = "2A12.20090929.67632.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

; 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 to plot.
height  =hdf4_file->heightLayerTop
;         printVarSummary(height)
;         printMinMax(height,True)
;         print(height)

CLUSTER          =hdf4_file->cluster(:,0,:,0) ; dim are ncluster,nlayer,nfindex,nspecies
;CLUSTER@_FillValue = -9999.9
     printVarSummary(CLUSTER)           
     printMinMax(CLUSTER,True)            

frzr=hdf4_file->freezingHeight
frzr@lon2d=longitude
frzr@lat2d=latitude
;     printVarSummary(frzr)
;     printMinMax(frzr,True)
frzr1            = short2flt(frzr)
frzr1@lon2d=longitude
frzr1@lat2d=latitude
frzr1@_FillValue = -9999
;                    printVarSummary(frzr1)
;                    printMinMax(frzr1,True)

surf_pptn=hdf4_file->surfacePrecipitation
surf_pptn@lon2d      = longitude
surf_pptn@lat2d      = latitude
surf_pptn@_FillValue = -9999.9
;     printVarSummary(surf_pptn)
;     printMinMax(surf_pptn,True)

conv_pptn=hdf4_file->convectPrecipitation
conv_pptn@lon2d      = longitude
conv_pptn@lat2d      = latitude
conv_pptn@_FillValue = -9999.9
;     printVarSummary(conv_pptn)
;     printMinMax(conv_pptn,True)

surf_rain   =hdf4_file->surfaceRain
surf_rain@lon2d      =longitude
surf_rain@lat2d      =latitude
surf_rain@_FillValue = -9999.9
;     printVarSummary(surf_rain)
;     printMinMax(surf_rain,True)

clud_wtrpath=hdf4_file->cloudWaterPath
clud_wtrpath@lon2d      =longitude
clud_wtrpath@lat2d      =latitude
clud_wtrpath@_FillValue = -9999.9
;     printVarSummary(clud_wtrpath)
;     printMinMax(clud_wtrpath,True)

clus_scale=hdf4_file->clusterScale(:,:,0)  
clus_scale@lon2d      =longitude
clus_scale@lat2d      =latitude
clus_scale@_FillValue = -9999.9
     printVarSummary(clus_scale)
     printMinMax(clus_scale,True)

profile_val_of_var0_at_layer0   = clus_scale*CLUSTER
                printVarSummary(profile_val_of_var0_at_layer0)

xwks=gsn_open_wks("png",file_name+".profile_val_of_var0_at_layer0")

res=True ; plot mods desired
res@cnFillOn=True ; enable contour fill
;res@gsnMaximize=True; make plot large
;res@gsnPaperOrientation =	"portrait" ; force portrait orientation
;res@cnLinesOn=False ; turn off contour lines
;res@cnLineLabelsOn       =  False; turn off contour line labels
;res@gsnSpreadColors=True ; use the entire color spectrum
;res@cnFillMode="RasterFill" ; faster
;res@lbOrientation="vertical" ; vertical labels
;res@cnMissingValFillPattern = 0 ; missing value pattern is set to "SolidFill"
;res@cnMissingValFillColor=0; white color for missing values
; Without the following, it will take a very long time to plot.
res@trGridType = "TriangularMesh" 
;res@cnLevelSelectionMode="ExplicitLevels";
;res@cnLevels=(/0, 0.1, 1.0, 10.0, 30.0/)
;res@cnLevels=(/0, 1.0, 10.0, 30.0/)
;gsn_define_colormap(xwks,"wh-bl-gr-ye-re") ; define colormap
gsn_define_colormap(xwks,"wh-bl-gr-ye-re") ; define colormap
;gsn_define_colormap(xwks,"sunshine_9lev") ; define colormap

res@tiMainString=file_name ; create title
;plot=gsn_csm_contour_map_ce(xwks,var0,res)
delete(plot) ; cleaning up resources used
delete(clus_scale)
delete(xwks)
delete(res)

delete(hdf4_file)
end

; References
; [1] http://disc.sci.gsfc.nasa.gov/precipitation/documentation/TRMM_README/TRMM_2A12_readme.shtml