; mls.ncl
;
; This file reads the remote MLS file via OPeNDAP and plots
; vertical profile (time x pressure).
; 
; This demo is based on the following example:
; http://www.ncl.ucar.edu/Applications/Scripts/eosdis_2.ncl
;
; and
;
; hdf5eos_3.ncl
; 
; Requires: NLC 5.1.1
;
; Author: Hyo-Kyung Lee (hyoklee@hdfgroup.org)
; 
; Copyright (C) 2009 The HDF Group
;
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/csm/contributed.ncl"

begin
url = "http://hdfdap.hdfgroup.uiuc.edu:8888/opendap/data/NASA/HDF5/HDF-EOS5/AURA/MLS/Swath/" 
fName = "MLS-Aura_L2GP-CO_v02-23-c02_2008d277.he5" 
if ( isfilepresent(url+fName) ) then
  f = addfile ( url + fName    , "r")
  vNames = getfilevarnames(f) 
end if

; Uncomment this to check if the remote file is opened successfully.
; print(f)
;
; Uncomment this to locate the L2gpValue.
; print(vNames)
;
; From the above, we know that the 12th variable is L2gpValue 
T = f->$vNames(12)$ 
; Re-scale the values.
T = T*1e7
T@units  = T@units + "*1e7"
lat = f->$vNames(8)$
lon = f->$vNames(5)$
; The below mapping assumes that the OPeNDAP hdf5_handler maps nTimes to time
; and Pressure to lev variable.
; Otherwise, you need to use numeric index.
lev = f->lev
time = f->time
ntime = dimsizes(time)

; This is useful to set manual contour levels. See "ManualLevels" below.
printVarSummary(T)
; The rest of the code is almost identical to hdf5eos_3.ncl
printMinMax(T, True)
print("==============")

; time units: NCL does not *yet* have a function analogous
; to ut_calendar for TAI93 (International Atomic Time). 
; Create an "elapsed time" variable.

    telapse           = (time - time(0))/60
    telapse@long_name = "Elapsed Time (minutes)"
    telapse@units     = "minutes since "+time(0)

;************************************************
; plot
;************************************************
pltType  = "ps"
pltName  = fName+".L2gpValue.ncl"

    wks   = gsn_open_wks (pltType,pltName)        ; open workstation
    gsn_define_colormap(wks,"amwg")               ; choose colormap
    dumb  = NhlNewColor(wks,0.7,0.7,0.7)          ; add gray to colormap
 
    res   = True
    res@gsnMaximize          = True               ; make ps/eps/pdf
    res@gsnSpreadColors      = True
    res@gsnSpreadColorEnd    = -2                 ; do not use gray for contours
   ;res@gsnPaperOrientation  = "portrait"         ; force portrait

    res@tiMainString         = fName

    res@cnFillOn             = True
    res@cnLinesOn            = False
    res@cnLineLabelsOn       = False
    res@cnFillMode           = "RasterFill"       ; faster
    res@cnRasterSmoothingOn  = True

    res@lbLabelAutoStride    = True
    res@lbOrientation        = "Vertical"

; Change these values based on min/max values of L2gpValue dataset. 
    res@cnLevelSelectionMode = "ManualLevels"     ; set manual contour levels
    res@cnMinLevelValF       = -6.               ; set min contour level
    res@cnMaxLevelValF       = 6.               ; set max contour level
    res@cnLevelSpacingF      =  1.               ; set contour spacing

    res@trYReverse           = True               ; reverse y-axis
    T&lev    = lev                ; assign "pressure" coordinates
    T&time     = telapse            ; assign temporal coordinates

    res@gsnYAxisIrregular2Linear = True
    res@tmYLFormat           = "f"                ; force minimal precision 
    res@tiYAxisString        = "Pressure (hPa)"

    res@gsnLeftString        = T@title

    plot  = gsn_csm_contour (wks,T(lev|:,time|:) , res)

;***************************************
; create trajectory plot
;***************************************

    mpres                    = True        ; Plot options desired.
    mpres@gsnFrame           = False       ; Don't advance the frame
    mpres@gsnMaximize        = True
    mpres@mpLandFillColor    = "gray70"    ; color of land
   ;mpres@gsnPaperOrientation= "portrait"  ; force portrait
    plot = gsn_csm_map_ce(wks,mpres)       ; Draw map   

    gsres                    = True        ; "Graphic Style" resources
    gsres@gsMarkerSizeF      = 10.0        ; Marker size
    gsres@gsMarkerThicknessF = 1.0         ; Marker thickness
    gsres@gsMarkerColor      = "Blue"      ; Marker color
    gsres@gsMarkerIndex      = 1           ; Marker style

    mpres@tiMainString       = fName

    gsn_polymarker(wks,plot,lon,lat,gsres) ; plot trajectory
    frame(wks)

end