;***********************************************
; plots the joint histograms of ISCCP, MODIS, and MISR
; following the NCL example: raster_6.ncl
; https://www.ncl.ucar.edu/Applications/Scripts/raster_6.ncl
;***********************************************

; Read in data and average. JH data are printed per output freq.
; There should be just as many files as the amount of processors used.

undef("modisP")
procedure modisP(f:list,res:logical)
local fq,cri,crl,im,iim,lim,var,wgt,liq,ice,x,y,icew,liqw
begin
  fq  = "modis_fq"
  crl = "modis_cotvreff_liq"
  cri = "modis_cotvreff_ice"
  
  var = dim_sum_n_Wrap(f[:]->joint_hist_fq,(/0/))
  printVarInfo(var,"var")
  wgt = dim_sum_n_Wrap(f[:]->joint_hist_fq_wgt,(/0/))
  printVarInfo(wgt,"wgt")
  ;Do final calculation of the weighted mean and then the monthly mean
  wgt@_FillValue = 0.0
  var@_FillValue = 0.0
  var = mask(var,wgt.eq.0.0,False)
  var = where(.not.ismissing(wgt),var/wgt,0.0)
  im = dim_avg_n_Wrap(var,(/0/))
  ; limit the decimal places, view data, and delete some atts
  im = decimalPlaces(im,2,True)
  im@units = "Cloud Fraction [%]"
  im@long_name = "Jan 2007 OpenIFS T255L91"
  ; 0 cannot be used as fillvalue when making contour
  delete_VarAtts(im,(/"standard_name","_FillValue"/))
  
  ; print out some useful info
  write_matrix (im,"7f5.2",False)
  printVarInfo(im,"im")
  
  liq = dim_sum_n_Wrap(f[:]->joint_hist_cotvreffliq,(/0/))
  printVarInfo(liq,"liq")
  liqw = dim_sum_n_Wrap(f[:]->joint_hist_cotvreffliq_wgt,(/0/))
  printVarInfo(liqw,"liqw")

  liqw@_FillValue = 0.0
  liq@_FillValue = 0.0

  liq = mask(liq,liqw.eq.0.0,False)
  liq = where(.not.ismissing(liqw),liq/liqw,0.0)

  lim = dim_avg_n_Wrap(liq,(/0/))                           
  lim = decimalPlaces(lim,2,True)
  lim@units = "Cloud Fraction [%]"
  lim@long_name = "Jan 2007 OpenIFS T255L91"
                                
  write_matrix (lim,"7f5.2",False)
  printVarInfo(lim,"lim")
  delete([/liq,liqw/])
  delete_VarAtts(lim,(/"standard_name","_FillValue"/))
      
  ice = dim_sum_n_Wrap(f[:]->joint_hist_cotvreffice,(/0/))
  printVarInfo(ice,"ice")
  icew = dim_sum_n_Wrap(f[:]->joint_hist_cotvreffice_wgt,(/0/))
  printVarInfo(icew,"icew")
  
  icew@_FillValue = 0.0
  ice@_FillValue = 0.0
  ice = mask(ice,icew.eq.0.0,False)
  ice = where(.not.ismissing(icew),ice/icew,0.0)
  iim = dim_avg_n_Wrap(ice,(/0/))
                            
  iim = decimalPlaces(iim,2,True)
  iim@units = "Cloud Fraction [%]"
  iim@long_name = "Jan 2007 OpenIFS T255L91"
  delete([/ice,icew/])
                     
  write_matrix (iim,"7f5.2",False)
  printVarInfo(iim,"iim")
  delete_VarAtts(iim,(/"standard_name","_FillValue"/))

  print("Setting up plot...")
  x = (/0.0,0.3,1.3,3.6,9.4,23.0,60.0,100.0/)
  res@sfXArray      = x
  res@tmXBMode      = "Explicit"
  res@tmXBValues    = x
  res@tmXBLabels    = (/"0.0","0.3","1.3","3.6","9.4","23.0","60.0"," "/)
  res@tiXAxisString = "Cloud Optical Depth"
  
  res@tmYLMode = "Explicit"
  res@trYReverse = True       
  y = (/1000.0,800.0,680.0,560.0,440.0,310.0,180.0,0.0/)
  res@sfYArray      = y
  res@tmYLValues    = y
  res@tmYLLabels    = (/"1000","800","680","560","440","310","180"," "/)
  res@tiYAxisString = "Cloud Top Pressure [hPa]"

  res@vpWidthF   = 0.9
  res@vpHeightF  = 0.6

  res@gsnYAxisIrregular2Linear = True
  
  wks = gsn_open_wks("png",fq) 
  print("Drawing plot modis fq...")
  plot = gsn_csm_contour(wks,im,res)

  wks = gsn_open_wks("png",cri) 
  print("Drawing plot modis cot v reff ice...")
  plot = gsn_csm_contour(wks,iim,res)

  wks = gsn_open_wks("png",crl) 
  print("Drawing plot modis cot v reff liq...")
  plot = gsn_csm_contour(wks,lim,res)
end

undef("isccpP")
procedure isccpP(f:list,res:logical)
local fq
begin
  fq = "isccp_fq"
  ; sum the files
  var = dim_sum_n_Wrap(f[:]->joint_hist_fq,(/0/))
  printVarInfo(var,"var")
  wgt = dim_sum_n_Wrap(f[:]->joint_hist_fq_wgt,(/0/))
  printVarInfo(wgt,"wgt")
  ; Do final calculation of the weighted mean and then the monthly mean
  wgt@_FillValue = 0.0
  var@_FillValue = 0.0
  var = mask(var,wgt.eq.0.0,False)
  var = where(.not.ismissing(wgt),var/wgt,0.0)
  im = dim_avg_n_Wrap(var,(/0/))
  ; limit the decimal places, view data, and delete some atts
  im = decimalPlaces(im,2,True)
  im@units = "Cloud Fraction [%]"
  im@long_name = "Jan 2007 OpenIFS T255L91"
  ; 0 cannot be used as fillvalue when making contour
  delete_VarAtts(im,(/"standard_name","_FillValue"/))
  
  ; print out some useful info
  write_matrix (im,"7f5.2",False)
  printVarInfo(im,"im")
  
  print("Setting up plot...")
  x = (/0.0,0.3,1.3,3.6,9.4,23.0,60.0,100.0/)
  res@sfXArray      = x
  res@tmXBMode      = "Explicit"
  res@tmXBValues    = x
  res@tmXBLabels    = (/"0.0","0.3","1.3","3.6","9.4","23.0","60.0"," "/)
  res@tiXAxisString = "Cloud Optical Depth"
  
  res@tmYLMode = "Explicit"
  res@trYReverse = True       
  y = (/1000.0,800.0,680.0,560.0,440.0,310.0,180.0,0.0/)
  res@sfYArray      = y
  res@tmYLValues    := y
  res@tmYLLabels    = (/"1000","800","680","560","440","310","180"," "/)
  res@tiYAxisString = "Cloud Top Pressure [hPa]"

  res@gsnYAxisIrregular2Linear = True
  res@vpWidthF   = 0.9
  res@vpHeightF  = 0.6
 
  wks = gsn_open_wks("png",fq) 
  print("Drawing plot...")
  plot = gsn_csm_contour(wks,im,res)
end

undef("misrP")
procedure misrP(f:list,res:logical)
local fq,var,wgt,im,x,y
begin
  fq = "misr_fq"
  ; sum the files
  var = dim_sum_n_Wrap(f[:]->joint_hist_fq,(/0/))
  printVarInfo(var,"var")
  wgt = dim_sum_n_Wrap(f[:]->joint_hist_fq_wgt,(/0/))
  printVarInfo(wgt,"wgt")
  
  wgt@_FillValue = 0.0
  var@_FillValue = 0.0
  var = mask(var,wgt.eq.0.0,False)
  var = where(.not.ismissing(wgt),var/wgt,0.0)
  im = dim_avg_n_Wrap(var,(/0/))
                                
  im = decimalPlaces(im,2,True)
  im@units = "Cloud Fraction [%]"
  im@long_name = "Jan 2007 OpenIFS T255L91"                           
  delete_VarAtts(im,(/"standard_name","_FillValue"/))
  
  write_matrix (im,"7f5.2",False)
  printVarInfo(im,"im")
  
  print("Setting up plot...")
  ; common to all cases
  x = (/0.0,0.3,1.3,3.6,9.4,23.0,60.0,100.0/)
  res@tmXBMode      = "Explicit"
  res@sfXArray      = x
  res@tmXBValues    := x
  res@tmXBLabels    = (/"0.0","0.3","1.3","3.6","9.4","23.0","60.0"," "/)
  res@tiXAxisString = "Cloud Optical Depth"

  y := (/-1.0,0.0,0.5,1.0,1.5,2.0,2.5,3.0,4.0,5.0,7.0,9.0,11.0,13.0,15.0,17.0,18.0/)
  res@tmYLMode      = "Explicit"
  res@sfYArray      := y
  res@tmYLValues    = y
  res@tmYLLabels    = (/"NR","0","0.5","1.0","1.5","2.0","2.5","3.0","4.0","5.0",\
  "7.0","9.0","11.0","13.0","15.0","17.0"," "/)
  res@tiYAxisString = "Cloud Top Height [km]"
  
  res@gsnYAxisIrregular2Linear = True

  res@vpWidthF   = 1.0
  res@vpHeightF  = 0.6
 
  wks = gsn_open_wks("png",fq) 
  print("Drawing plot misr ctp vs hgt...")
  plot = gsn_csm_contour(wks,im,res)  
end

undef("lidarP")
procedure lidarP(f:list,res:logical)
local sr
begin
  sr = "lidar_sr_vs_hgt"
  var = dim_sum_n_Wrap(f[:]->sr_vs_hgt,(/0/))
  printVarInfo(var,"var")
  wgt = dim_sum_n_Wrap(f[:]->sr_vs_hgt_wgt,(/0/))
  printVarInfo(wgt,"wgt")
  
  wgt@_FillValue = 0.0
  var@_FillValue = 0.0
  var = mask(var,wgt.eq.0.0,False)
  var = where(.not.ismissing(wgt),var/wgt,0.0)
  im = dim_avg_n_Wrap(var,(/0/))
                                
  im = decimalPlaces(im,2,True)
  im@units = "Cloud Fraction [%]"
  im@long_name = "Jan 2007 OpenIFS T255L91"
                                
  write_matrix (im,"7f5.2",False)
  printVarInfo(im,"im")
  delete([/var,wgt/])
  delete_VarAtts(im,(/"standard_name","_FillValue"/))
  im&sr = decimalPlaces(im&sr,2,True)
  im&hgt = decimalPlaces(im&hgt * 1.0e-3,2,True)
      
  print("Setting up plot...")

  x = (/0.0,0.01,1.2,3.0,5.0,7.0,10.0,15.0,20.0,25.0,30.0,40.0,50.0,60.0,80.0,81.0/)
  res@sfXArray      = x
  res@tmXBMode      = "Explicit"
  res@tmXBValues    = x
  res@tmXBLabels    = (/"0.0","0.01","1.2","3.0","5.0","7.0","10.0","15.0","20.0",\
                        "25.0","30.0","40.0","50.0","60.0","80.0"," "/)
  res@tiXAxisString = "Back-scattering Ratio"
  
  ;res@tmYLMode      = "Explicit"
  ;res@trYReverse = True
  ;y = (/240,720,1200,1680,2160,2640,3120,3600,4080,4560,5040, \
  ;      5520,6000,6480,6960,7440,7920,8400,8880,9360,9840,10320, \
  ;      10800,11280,11760,12240,12720,13200,13680,14160,14640,15120, \
  ;      15600,16080,16560,17040,17520,18000,18480,18960/) * 1.0e-3
  ;y = decimalPlaces(y,1,True)
  ;res@sfYArray      := y
  ;res@tmYLValues    = y
  ;res@tmYLLabels    = (/"0.2","0.7","1.2","1.6","2.2","2.6","3.1",\
  ;                      "3.6","4.1","4.6","5.0","5.5","6.0","6.5",\
  ;                      "7.0","7.4","7.9","8.4","8.9","9.4","9.8","10.0",\
  ;                      "10.8","11.3","11.8","12.2","12.7","13.2","13.7",\
  ;                      "14.1","14.6","15.1","15.6","16.1","16.6","17.0",\
  ;                      "17.5","18.0","18.5","19.0"/)    
  ; Control the figure shape
  res@vpWidthF   = 0.9
  res@vpHeightF  = 0.7
  res@tiYAxisString = "Height [km]"
  
  wks = gsn_open_wks("png",sr) 
  print("Drawing plot...")
  plot = gsn_csm_contour(wks,im,res)
end

; ----
; Main
; ----
begin
  idir="/Users/marston/Desktop/hist_data/"
  isccp = idir+"isccp_*_fq_jhist.nc"
  modis = idir+"modis_*_fq_jhist.nc"
  misr  = idir+"misr_*_fq_jhist.nc"
  lidar = idir+"lidar_*_sr_vs_hgt.nc"

  cases = [/isccp,modis,misr,lidar/]

  do case = 0,3

    res                    = True         
    res@cnFillOn           = True        
    res@cnFillMode         = "RasterFill" 
    res@cnLinesOn          = False        
    res@cnLineLabelsOn     = False        
    res@cnFillPalette      = "BlAqGrYeOrReVi200" ;MPL_viridis"
    res@lbLabelBarOn       = True
    res@lbOrientation      = "vertical"    
    res@cnLabelBarEndStyle = "IncludeMinMaxLabels"
    res@lbBoxEndCapStyle   = "TriangleBothEnds"
    res@pmLabelBarOrthogonalPosF = -0.02
                                ; Control the figure shape
    res@vpWidthF   = 0.9
    res@vpHeightF  = 0.5
                                ; remove ticks not needed  
    res@tmYROn = False
    res@tmXTOn = False
    res@cnMonoFillPattern = True
    res@tmXBLabelFontHeightF = 0.014
    res@tmYLLabelFontHeightF = 0.014
    
    res@gsnPaperOrientation = "auto"
    res@gsnMaximize         = True
    
    cas := cases[case]
    print(" ")
    ; read files
    files = systemfunc ("ls "+cas)
    print("Found "+dimsizes(files)+" files!")
    print("Opening files in joined list...")
    nfiles = dimsizes(files)
    ; addfiles to list and join it
    f := addfiles(files,"r")
    ListSetType (f, "join")
    if(case.eq.0) then
      print("Processing and plotting data for isccp")
      isccpP(f,res)
    end if
    if(case.eq.1) then
      print("Processing and plotting data for modis")
      modisP(f,res)
    end if
    if(case.eq.2) then
      print("Processing and plotting data for misr")
      misrP(f,res)
    end if
    if(case.eq.3) then
      print("Processing and plotting data for lidar")
      lidarP(f,res)
    end if
    print("Complete!")
    delete(res)
  end do
end