;;dimensions:
;;	Scanline = 116 ;
;;	Field_of_view = 30 ;
;;	Channel = 20 ;
;;	P_Layer = 100 ;
;;	P_Level = 101 ;
;;	Qc_dim = 4 ;

;;      float Latitude(Scanline, Field_of_view) ;
;;              Latitude:long_name = "Latitude of the view (-90,90)" ;
;;      float Longitude(Scanline, Field_of_view) ;
;;              Longitude:long_name = "Longitude of the view (-180,180)" ;
;;      float PVapor(Scanline, Field_of_view, P_Layer) ;
;;              PVapor:long_name = "Water vapor profile in g/kg" ;
;;              PVapor:scale = 1. ;
;;      float Player(P_Layer) ;
;;              Player:description = "Pressure for each layer in mb" ;
;;      float Plevel(P_Level) ;
;;              Plevel:description = "Pressure for each level in mb" ;
;;---
;;      double ScanTime_UTC(Scanline) ;
;;              ScanTime_UTC:long_name = "Number of seconds since 00:00:00 UTC" ;
;;      short ScanTime_dom(Scanline) ;
;;              ScanTime_dom:long_name = "Calendar day of the month 1-31" ;
;;      short ScanTime_doy(Scanline) ;
;;              ScanTime_doy:long_name = "julian day 1-366" ;
;;      short ScanTime_hour(Scanline) ;
;;              ScanTime_hour:long_name = "hour of the day 0-23" ;
;;      short ScanTime_minute(Scanline) ;
;;              ScanTime_minute:long_name = "minute of the hour 0-59" ;
;;      short ScanTime_month(Scanline) ;
;;              ScanTime_month:long_name = "Calendar month 1-12" ;
;;      short ScanTime_second(Scanline) ;
;;              ScanTime_second:long_name = "second of the minute 0-59" ;
;;      short ScanTime_year(Scanline) ;
;;              ScanTime_year:long_name = "Calendar Year 20XX" ;
;;---
;;   // global attributes:
;;              :missing_value = -999 ;
;;              :cdf_version = 4. ;

      minlev  = 0
      maxlev  = 100   
      varName = "PVapor"

      pltType = "png"
      pltName = "MIRS."+varName
      pltDir  = "./"

    ;-- define file name
    ;;fili  = "NPR.MIRS.V7.SND.AAMH.NP.D18289.S0854.E1047.B4992829.NS.nc"
      fili  = "NPR.MIRS.V7.SND.AAMH.NP.D17165.S2035.E2230.B4303234.NS.nc"
    
    ;-- open file and read variables 
      f     =   addfile(fili, "r")  
      Player=   f->Player                ; Player(P_Layer)    ;-- [*]
      iPlay =   ind(Player.ge.minlev .and. Player.le.maxlev)
                                         ;--  Scanline ==> Time ; Read only desired levels
      var   =   f->$varName$(:,:,iPlay)  ;-- [Scanline | 116] x [Field_of_view | 30] x [P_Layer | ...]
      var@_FillValue = -999.0            ;-- specified by the global attributes       
printVarSummary(var)
printMinMax(var,0)
print("---")

      lat2d = f->Latitude                ; [Scanline | 116] x [Field_of_view | 30]
printVarSummary(lat2d)
printMinMax(lat2d,0)          
print("---")
      lon2d = f->Longitude
printVarSummary(lon2d)
printMinMax(lon2d,0)       
print("---")
                                         ; ScanTime ==> Scanline
      dom   = f->ScanTime_dom            ; day of month  
      doy   = f->ScanTime_doy            ; julian day 1-366                     
      hour  = f->ScanTime_hour
      minute= f->ScanTime_minute
      month = f->ScanTime_month
      sec   = f->ScanTime_second
      year  = f->ScanTime_year

      ymdhm = year*100000000l + 1000000l*month + 10000l*dom + 100l*hour +minute
      YMDHM = year+"/"+sprinti("%0.2i",month)+"/"+sprinti("%0.2i", dom) \
                  +"/"+sprinti("%0.2i",hour) +"/"+sprinti("%0.2i", minute)
print(ymdhm+"  "+YMDHM)  
print("---")
     
      minlat = -87 ; or some val
      maxlat =  76 ; or some val

      dimvar = dimsizes(var)
      nScan  = dimvar(0)
      nField = dimvar(1)
      nLayer = dimvar(2)
     
      pLay   =   Player(iPlay)
      pLay@long_name = "Layer Pressure"
      pLay@units = "hPa"
      pLay!0     = "pLay"
      pLay&pLay  =  pLay
printVarSummary(pLay)
printMinMax(pLay,0)
print("---")

      var!2      = "pLay"
      var&pLay   =  pLay
printVarSummary(var)
printMinMax(var,0)
print("---")
                              ; Use to set global contour min/max
      opt_stat = True
      opt_stat@PrintStat = True
      stat_var = stat_dispersion(var, opt_stat )  
                                                  
      maxvlev = 20            ; use nice_mnmxintvl to get 'nice' values
     ;mnmxint = nice_mnmxintvl( stat_var(3), stat_var(13), maxvlev, False) ; low/high dectiles
      mnmxint = nice_mnmxintvl( stat_var(3), stat_var(14), maxvlev, False) ; low dectile & max value
    
    ;-- set resources
      res                          =  True
    
      res@cnFillPalette            = "BlAqGrYeOrReVi200"
      res@cnFillOn                 =  True       ;-- turn on color fill
      res@cnFillMode               = "RasterFill"       ; Raster Mode
      res@cnLinesOn                =  False      ;-- Turn off contour lines
      res@cnLineLabelsOn           =  False      ;-- turns off contour line labels
      res@cnInfoLabelOn            =  False      ;-- turns off contour info label
      res@tiYAxisString            =  pLay@long_name + " [hPa] "
      res@trYReverse               =  True       ;-- reverses y-axis
      res@tmYLFormat               =  "f"        ; No unnecessary zeroes
      res@lbOrientation            = "vertical"  ;-- vertical label bar
      res@lbLabelFontHeightF       =  0.0150     ; change label size
 
      res@cnLevelSelectionMode     = "ManualLevels"   ; common contour levels for all plots
      res@cnMinLevelValF           = mnmxint(0)
      res@cnMaxLevelValF           = mnmxint(1)
      res@cnLevelSpacingF          = mnmxint(2)
    
    ;-- lat/lon label resources

      nLabels                      =  9
      res@tmXBValues               = fspan(0,nField-1,nLabels) 
      res@tmXBLabels               = new(nLabels,"string")
      res@tmXBMode                 = "Explicit"
      res@tmXBLabelFontHeightF     = 0.0150

      do ns=0,0   ; nScan-1

   ;-- define workstation
         pltPath = pltDir + pltName+"_"+sprinti("%0.3i",ns)

         wks = gsn_open_wks(pltType,pltPath)
         res@tiMainString          = "MIRS: Scan="+sprinti("%0.3i",ns)+": "+YMDHM(ns)

         do i=0,nLabels-1                         ;-- generate the lat/lon labels
            xx = lon2d(ns,i)
            yy = lat2d(ns,i)
            res@tmXBLabels(i) = sprintf("%5.1f", yy)+"~C~"+sprintf("%5.1f", xx)
         end do

         varPlt = var(Scanline|ns,pLay|:,Field_of_view|:)  ; clarity; reorder to (y,x)
         plot   = gsn_csm_contour(wks,varPlt,res)
      end do