;***********************************************
; 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.
begin
  idir="/Users/marston/Desktop/fq_data/"
  isccp = idir+"isccp_*_fq_jhist.nc"
  modis = idir+"modis_*_fq_jhist.nc"
  misr = idir+"misr_*_fq_jhist.nc"

  fnames = (/"isccp_fq","modis_fq","misr_fq"/)
  mylist = [/isccp,modis,misr/]

  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
  res@trYReverse = True       
  ; Control the figure shape
  res@vpWidthF   = 0.8
  res@vpHeightF  = 0.5
  ; remove ticks not needed  
  res@tmYROn = False
  res@tmXTOn = False
 
  res@gsnPaperOrientation = "auto"
  res@gsnMaximize         = True   
  
  do case = 0,2
    cas = mylist[case]
    print(" ")
    print("Processing and plotting data for = "+cas)
    ; 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")

    ; 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 [%]"

    ; print out some useful info
    write_matrix (im,"7f5.2",False)
    printVarInfo(im,"im")

    ; 0 cannot be used as fillvalue when making contour
    delete_VarAtts(im,(/"standard_name","_FillValue"/))
    im&tau = decimalPlaces(im&tau,1,True)

    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@sfXArray      := x
    res@tmXBMode      = "Explicit"
    res@tmXBValues    = x
    res@tmXBLabels    = ""+x
    res@tiXAxisString = "Cloud Optical Depth"
    
    if(case.lt.2) then ; modis and isccp
      y = (/1000.0,800.0,680.0,560.0,440.0,310.0,180.0,0.0/)
      im&pres = decimalPlaces(im&pres,0,True)
      res@sfYArray      := y
      res@tmYLMode      = "Explicit"
      res@tmYLValues    = y
      res@tmYLLabels    = (/"1000","800","680","560","440","310","180"," "/)
      res@tiYAxisString = "Cloud Top Pressure [hPa]"
    else
      ; misr
      y := (/-99.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,99.0/)
      res@sfYArray      := y
      res@tmYLMode      = "Explicit"
      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]"      
    end if

    res@cnMonoFillPattern = True
    ;res@tmXBMajorLengthF     = 0.0
    ;res@tmYLMajorLengthF     = 0.0
    ;res@tmXBMinorOn          = False
    ;res@tmYLMinorOn          = False
  
    delete([/files,f,var,wgt/])
    wks = gsn_open_wks("png",fnames(case)) 
    
    print("Drawing plot...")
    plot = gsn_csm_contour(wks,im,res)     
    delete([/wks,plot,im/])
    print("Complete!")
  end do
end