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

file_name = "CAL_LID_L2_05kmALay-Prov-V3-40.2019-10-09T00-12-56ZD.hdf"
hdf4_file = addfile("/global/noscrub/Partha.Bhattacharjee/NRT_NGACv2/GSD_Chem/Calipso/20191009/"+file_name, "r") 

fig_diri = "/global/noscrub/Partha.Bhattacharjee/NRT_NGACv2/GSD_Chem/Calipso/" ;directory that image gets saved in
fig_save = fig_diri + "CALIPSO_Backscatter_Test_2019-19-09T00-12-56ZD"; will save this image in the directory specified (fig_diri)
wks = gsn_open_wks("png", fig_save) ; open workstation & send graphics to PNG file. Saves graphics with file as indicated above. 

ext = hdf4_file->Column_Optical_Depth_Aerosols_532
printVarSummary(ext) ; Info:

  ext!0 = "first"
  ext!1 = "second"
  ext1 = ext(:,::-1)
  ext2  = ext1(second|:,first|:)
  ext2@_FillValue = -9999
  ext2@units      = "km~S~-1~NN~ sr~S~-1~NN~"

   print("ext2 (2018):  min="+min(ext2)+"  max="+max(ext2))
  printVarSummary(ext2)

Latitude = hdf4_file->Latitude
printVarSummary(Latitude)
Longitude = hdf4_file->Longitude
printVarSummary(Longitude)
xlabel = sprintf("%.2f",Latitude)+"~C~"+sprintf("%.2f",Longitude) ; format labeling on x-axis

   time    = hdf4_file->Profile_Time
   time@units = "seconds since 1993-01-01 00:00"
   tstring = ut_string(time(:,0),"%Y-%N-%D %H:%M:%S")
   
;----------------------------------------------
;  GET ALTITUDE VALUES; Figure out a way to extract metadata
;----------------------------------------------

   hgt    = asciiread("./altitudes.txt",-1,"float") ; hgt = height
 
;;    hgt = hdf4_file->Midlayer_Pressure
  
   printVarSummary(hgt) ;Info:

    hgt    = hgt(::-1) ; Reverse this array 
     hgt!0  = "hgt" ; first dimension being referenced, which is altitude.
     hgt@long_name = "Altitude, km"
     hgt@units     = "km"
     printVarSummary(hgt)


    ext3d = onedtond(ndtooned(ext2),(/3728,3728,399/))
    print(ext3d(360:361,698:699,40:41))
    ext3d!0 = "lat"
   ext3d!1 = "lon"
   ext3d!2 = "hgt"
  lon    = fspan(0,dimsizes(Longitude(:,0)) - 1,1)
  lat    = fspan(0,dimsizes(Latitude(:,0)) - 1,1) 
  
  lat@longname = "latitude"
  lat@units = "degrees_north"
  lon@longname ="longitude"
  lon@units = "degrees_east"

  print("ext3d (2018):  min="+min(ext3d)+"  max="+max(ext3d))
  printVarSummary(ext3d)


;---------------------------------------------- 
; PLOTTING RESOURCES
;---------------------------------------------- 

   setvalues NhlGetWorkspaceObjectId()
        "wsMaximumSize": 300000000
   end setvalues

   colorMap = read_colormap_file("BlAqGrYeOrRe")
   colorMap(0,:) = (/ 0., 0., 0., 1./)   
   res = True
   res@cnFillOn            = True            ; color plot desired
   res@cnLinesOn           = False
   res@cnFillPalette 	   = colorMap
   res@cnRasterModeOn      = True
   res@gsnAddCyclic        = False
   res@gsnMaximize         = True
   res@gsnLeftString       = "UTC: "+tstring(0)+" to "+tstring(dimsizes(tstring)-1)
   res@tiMainString        = "Attenuated Depolarization Ratio"
   res@tiMainFontHeightF   = 0.025
   res@tiYAxisFontHeightF  = 0.02
   res@sfYArray = hgt
  
;---------------------------------------------
;  HERE, SET SPECIAL LABEL HEIGHTS AND ORIENTATION
;---------------------------------------------
   res@lbLabelFontHeightF = 0.01
   res@lbLabelAngleF = 0.0
   res@lbOrientation = "vertical"
   res@vpWidthF = 1.0
   
;---------------------------------------------
;  HERE YOU CAN SET A MAXIMUM ALTITUDE
;---------------------------------------------

   res@trYMaxF = 10  ; THIS IS IN KILOMETERS. Needed? Check once plot has been constructed
   
;---------------------------------------------
;  SET SPECIAL LEVELS
;---------------------------------------------

  res@cnLevelSelectionMode = "ExplicitLevels" 
  res@cnLevels             = (/ 0.0001, 0.0002, 0.0003, 0.0004, 0.0005, 0.0006, 0.0007, 0.0008, 0.0009, \
                                  0.001, 0.0015,  0.002, 0.0025,  0.003, 0.0035,  0.004, 0.0045,  0.005, \
                                 0.0055,  0.006, 0.0065,  0.007,  0.0075, 0.008, \ 
                                   0.01,   0.02,   0.03,   0.04,   0.05,   0.06,   0.07,   0.08,   0.09,  0.1 /) 
                                   
  ;levels from example NCL script. Appears to be correct from Browse Images section in Subsetter Tool                                

;---------------------------------------------
;  HERE IS WHERE YOU STRIDE THE DATA AND PLOT LATITUDE AND LONGITUDE ON X-AXIS
;---------------------------------------------

  ;tm = tickmark
   res@tmXTOn = False ; turns top tick marks off
   res@tmXBMode = "Explicit" ;uses the array resource tmXBValues to determine the coordinate positions of the tick marks
   res@tmXBValues = xcoord(::stride)  ;(::stride) is used to subsample an array. subsamples xcoord. 
   res@tmXBLabels = xlabel(::stride,0)
   res@tmXBLabelFontHeightF = 0.01 ; toggle this if you want
   
;---------------------------------------------
;  HERE, SET SPECIAL LABEL HEIGHTS AND ORIENTATION
;---------------------------------------------

   res@lbLabelFontHeightF = 0.01
   res@lbLabelAngleF = 0.0
   res@lbOrientation = "vertical"
   res@vpWidthF = 1.0
   
;---------------------------------------------
;  PLOT
;---------------------------------------------
  
   res@gsnDraw = True
   res@gsnFrame = False
   plot = gsn_csm_contour(wks,ext3,res)
   frame(wks)

;delete([/Latitude,Longitude, ext, xcoord, xlabel/]) -- probably important once loop is constructed





end