; The CALIPSO Level 2 lidar Vertical Feature Mask (VFM) data product describes the vertical and 
; horizontal distribution of cloud and aerosol layers observed by the CALIPSO lidar.
;=====================================================================
; https://www-calipso.larc.nasa.gov/resources/calipso_users_guide/data_summaries/vfm/ 
;=====================================================================
  dirCAL = "/Users/shea/Data/CALIPSO/CAL_LID/"
  filCAL = systemfunc("cd "+dirCAL+" ; ls CAL_LID_L2_VFM-ValStage1-V3-01.2010*hdf")
  nfil   = dimsizes(filCAL)
  print(filCAL)
  print("=============")

  nskip   = 16    ; two bytes; bit size of types short and ushort

  do nf=0,nfil-1  ; loop over each file
     print("--------------------------------------------------------------------")
     print("File Name: "+filCAL(nf))
     print("--------------------------------------------------------------------")

;---Extract variables on the file
;---NCL's reassignment operator [ := ] allows for different variable sizes for each file

     f       = addfile(dirCAL+filCAL(nf),"r")

     time   := f->Profile_Time(:,0)               ; (time)
     time!0  = "time"
     time@long_name = time@hdf_name               ; units are terrible: "yymmdd.frac"
     ntim    = dimsizes(time)

     lat    := f->Latitude(:,0)                   ; (time)                       
     lat!0   = "time"
     lat@long_name = lat@hdf_name 
     lat@units     = "degrees_north"              ; units recognized by NCL

     lon    := f->Longitude(:,0)                  ; (time)             
     lon!0   = "time"
     lon@long_name = "Longutude"
     lon@units     = "degrees_east"               ; units recognized by NCL

     dn     := f->Day_Night_Flag(:,0)             ; (time)
     dn!0    = "time"
     dn@long_name = "Day_Night_Flag"
     dn@units     = "0=day, 1=night"

     lwm    := f->Land_Water_Mask(:,0)            ; (time)
     lwm!0   = "time"
     lwm@long_name = lwm@hdf_name      
     lwm@units    := (/ "0=shallow ocean"       \
                      , "1=land"                \
                      , "2=coastlines"          \
                      , "3=shallow inland water"\
                      , "4=intermittent water"  \
                      , "5=deep inland water"   \
                      , "6=continental ocean"   \
                      , "7=deep ocean"          /)


     sp     := f->Spacecraft_Position             ; (time,3) 
     sp!0    = "time"
     sp!1    = "sat"
     sp@long_name = sp@hdf_name
     sp@_FillValue= sp@fillvalue                  ; NCL requires _FillValue

     fcf    := f->Feature_Classification_Flags    ; (time,layer)  ; ushort ; ushort (16 bits <==> 2 bytes)
     fcf@long_name = fcf@hdf_name
     fcf@info      = "cloud vs. aerosol vs. stratospheric layer" 
     fcf!0  = "time"
     fcf!1  = "layer"                             ; 5515
     
     dimfcf  = dimsizes(fcf)
     nlay    = dimfcf(1)                          ; 5515

;---Unpack 'fcf' bits
;---Details: 
;https://www-calipso.larc.nasa.gov/resources/calipso_users_guide/data_summaries/vfm/#feature_classification_flags

     ftyp   := dim_gbits (fcf, 0, 3, nskip, nlay)   ; Bits 1-3; Feature Type; 
     ftyp!0  = "time"
     ftyp!1  = "layer"                             
     ftyp@long_name = "Feature Type"    
     ftyp@units = (/ "0=invalid (bad or missing data)"  \
                   , "1=clear air"                      \
                   , "2=cloud"                          \    ; see ctyp below
                   , "3=aerosol"                        \    ; see atyp below
                   , "4=stratospheric feature"          \
                   , "5=surface"                        \
                   , "6=subsurface"                     \
                   , "7=no signal (totally attenuated)" /)
printVarSummary(ftyp)
print("ftyp: min="+min(ftyp)+"  max="+max(ftyp))
print("------------")

     fqul   := dim_gbits (fcf, 3, 2, nskip, nlay)   ; Bits 4-5; Feature Type Quality 
     fqul@long_name =  "Feature Type Quality"
     fqul!0  = "time"
     fqul!1  = "layer"                             

     iwp    := dim_gbits (fcf, 5, 2, nskip, nlay)  ; Bits 6-7; Ice/Water Phase
     iwp@long_name = "Ice/Water Phase"
     iwp!0   = "time"
     iwp!1   = "layer"                             

     iwpq   := dim_gbits (fcf, 7, 2, nskip, nlay)  ; Bits 8-9; Ice/Water Phase Quality
     iwpq@long_name = "Ice/Water Phase Quality"
     iwpq!0  = "time"
     iwpq!1  = "layer"                             

     styp   := dim_gbits (fcf, 9, 3, nskip ,nlay)  ; sub-type ==> Bits 10-12 ==> ushort     
     styp@long_name = "Feature Sub-Type"
     styp!0  = "time"
     styp!1  = "layer"                             
printVarSummary(styp)
print("styp: min="+min(styp)+"  max="+max(styp))

     caqul  := dim_gbits (fcf,12, 1, nskip, nlay)  ; Bit 13; Cloud / Aerosol /PSC Type QA
     caqul@long_name = "CLoud/Aerosol/PSC Phase Quality"
     caqul!0 = "time"
     caqul!1 = "layer"                             

     havg   := dim_gbits (fcf,13, 3, nskip, nlay) ; Bits 14-16; Horizanal Averaging
     havg@long_name = "Horizanal Averaging"
     havg!0  = "time"
     havg!1  = "layer"                             

;---Create explicit cloud and aerosol arrays

     fill  = totype(99,typeof(styp))
     CTYP := where(ftyp.eq.2, styp, fill)                    ; cloud type
     CTYP@_FillValue = fill

     ctyp := toshort(CTYP)                                   ; prefer "short"
     ctyp@long_name = "Cloud Type"
     ctyp!0  = "time"
     ctyp!1  = "layer"                             
     ctyp@units = (/ "1=low overcast, opaque"     \
                   , "2=transition stratocumulus" \
                   , "3=low, broken cumulus"      \
                   , "4=altocumulus (transparent)"\
                   , "5=altostratus (opaque)"     \
                   , "6=cirrus (transparent)"     \
                   , "7=deep convective (opaque)" /)

     ATYP := where(ftyp.eq.3, styp, fill)                    ; aerosol type
     ATYP@_FillValue = fill

     atyp := toshort(ATYP)                                   ; prefer "short"
     atyp@long_name = "Aerosol Type"
     atyp!0  = "time"
     atyp!1  = "layer"                             
     atyp@units = (/ "0=not determined"      \
                   , "1=clean marine"        \
                   , "2=dust"                \
                   , "3=polluted continental"\
                   , "4=clean continental"   \
                   , "5=polluted dust"       \
                   , "6=smoke"               \
                   , "7=other"               /)


printVarSummary(ctyp)
print("ctyp: min="+min(ctyp)+"  max="+max(ctyp))
print("-----")

printVarSummary(atyp)
print("atyp: min="+min(atyp)+"  max="+max(atyp))
print("-----")

 
    ;ad    := where(lat(1:).gt.lat(0:nlat-2), 1, -1) ; ascending/descending     
    ;ad!0   = "time"
    ;nad    = dimsizes(ad)
    ;npass  = num(ad(0:nad-2).ne.ad(1:)) + 1
  end do