; Concepts illustrated:
;   - Reading TRMM netCDF file                          
;   - Read a variable of type short and convert to float          
;   - Use NCL's distributed 1x1 land-sea mask top create a mask file 
;     CAVEAT: landsea_mask uses a 1x1 basemao. TRMM 3B42 is 0.25 
;   - Writing mask array to a NetCDF
;   - Use the mask to plot over land and lakes only  
;***************************************************************
  netCDF   = False
  PLOT     = True

;--- landsea.nc is distributed with NCL
;--- CAVEAT: landsea.nc has 1.0  degree resolution
;--- TRMM:   0.25 degree resolution
 
;--- Read TRMM grid  

  f = addfile("yearsum2000.nc","r")
  TRMM = f->pre

  TRMM_reorder = TRMM(time|:,latitude|:,longitude|:) 
  ;printVarSummary(TRMM_reorder)

if (any(isnan_ieee(TRMM_reorder))) then
      value = 1.e20
      replace_ieeenan (TRMM_reorder, value, 0)
      TRMM_reorder@_FillValue = value
end if
 
  ;printMinMax(TRMM_reorder,0)
  ;print("---")

;--- landsea.nc is distributed with NCL
;--- read 1x1 land sea mask basemap file
;--- create mask array; NOTE basemap is 1 deg while TRMM is 0.25 deg                      
;--- make mask contain byte values for smaller file

  f_lsm    = addfile("$NCARG_ROOT/lib/ncarg/data/cdf/landsea.nc","r")  

  lsmTRMM  = tobyte( landsea_mask(f_lsm->LSMASK,TRMM_reorder&latitude,TRMM_reorder&longitude)   )
  lsmTRMM@long_name  = "Land-Sea Mask: TRMM 3B42"
  ;printVarSummary(lsmTRMM)
  copy_VarCoords(TRMM_reorder(0,:,:), lsmTRMM)
  ;printMinMax(lsmTRMM,0)
  ;print("---")

;====================== netCDF ==============================================
if (netCDF) then
  dirNC = "./"
  filNC = "TRMM_3B42.LandSeaMask.nc"
  pthNC = dirNC + filNC
  system("/bin/rm -f "+pthNC)       ; remove any pre-existing file
  ncdf = addfile(pthNC, "c")        ; open output netCDF file

  fAtt               = True            ; assign file attributes
  fAtt@title         = "NCL 1-degree basemap onto 0.25 deg TRMM 3B42 grid"
  fAtt@Conventions   = "None"
  fAtt@creation_date = systemfunc ("date")
  fileattdef( ncdf, fAtt )            ; copy file attributes, Defines global (file) attributes associated with a supported file.
  ncdf->LSMASK_3B42  = lsmTRMM
end if

;======================= PLOT ===============================================
if (PLOT) then
  wks = gsn_open_wks("png","3B42_mask")  ; send graphics to PNG file
  
setvalues NhlGetWorkspaceObjectId() 
    "wsMaximumSize" : 300000000
  end setvalues 

  res = True
  res@mpFillOn         = False                ; do not color-fill the map
  res@mpPerimDrawOrder = "PostDraw"           ; draw the map perimeter last
  res@mpMaxLatF        = max(TRMM_reorder&latitude)
  res@mpMinLatF        = min(TRMM_reorder&latitude)

  res@cnFillOn         = True                 ; turn on color fill
  res@cnLinesOn        = False                ; turn off the contour lines
  res@cnLineLabelsOn   = False                ; turn off the contour line labels
  res@cnFillMode       = "RasterFill"         ; use raster fill 
  res@cnRasterSmoothingOn = True        
  res@cnFillPalette    = "wh-bl-gr-ye-re"
  res@cnLevelSelectionMode = "ExplicitLevels" ; explicitly set the levels via cnLevels
  res@cnLevels         = (/1.,2.,3.,4./)      ; set the levels
  res@cnFillColors     = (/60,100,20,140,5/)  ; set the colors that will be used to color fill

  res@lbLabelStrings   = ispan(0,4,1)         ; labels for the labelbar boxes
  res@lbLabelAlignment = "BoxCenters"         ; put the labels in the center of the label bar boxes
  res@lbTitleString    = "0=ocean, 1=land, 2=lake, 3=small island, 4=ice shelf"  ; labelbar title
  res@lbTitleFontHeightF       = 0.0125       ; labelbar title font height
  res@pmLabelBarOrthogonalPosF = .125         ; move whole thing down
  
  res@gsnRightString = "CAVEAT: 1 deg basemap used: TRMM 0.25 deg"
  plt  = gsn_csm_contour_map(wks,lsmTRMM,res)

  plot = new(3,graphic)
  res2 = True
  res2@gsnDraw          = False                ; do not draw the plot
  res2@gsnFrame         = False                ; do not advance the frame

  res2@mpFillOn         = False                ; do not color-fill the map
  res2@mpPerimOn        = True                 ; turn the map perimeter on
  res2@mpPerimDrawOrder = "PostDraw"           ; draw the map perimeter last
  res2@mpMaxLatF        = max(TRMM_reorder&latitude)
  res2@mpMinLatF        = min(TRMM_reorder&latitude)

  res2@cnLinesOn        = False                ; turn off the contour lines
  res2@cnLineLabelsOn   = False                ; turn off the contour line labels
  res2@cnFillOn         = True                 ; turn on color fill
  res2@cnLevelSelectionMode = "ExplicitLevels"              
  res2@cnLevels         = (/0.1,1.0,1.5,2.0,2.5,5,7.5,10,15,20/) ; "mm/day" 
  res2@cnLevels         = 365*res2@cnLevels    ; annual total
  res2@cnFillPalette    = (/"Snow","PaleTurquoise","PaleGreen" \
                           ,"SeaGreen3" ,"Yellow", "Orange","HotPink" \ ; contour colors
                           ,"Red","Violet", "Purple", "Brown"  /)       ; one more color than contour levels
  res2@cnMissingValFillPattern = 0             ; make 'missing' black
  res2@cnMissingValFillColor   = "black"       ; "background"

  res2@lbLabelBarOn      = False           ; turn off individual cb's

  nt = 0  ; 1st timestep
  plot(0) = gsn_csm_contour_map(wks,TRMM_reorder(nt,:,:),res2)
  
; apply the Land-Sea-Mask 
  prcMASK = where(lsmTRMM.eq.1 .or. lsmTRMM.eq.2, TRMM_reorder(nt,:,:), TRMM_reorder@_FillValue) ; Land only
  copy_VarMeta(TRMM_reorder(nt,:,:), prcMASK)
  plot(1) = gsn_csm_contour_map(wks,prcMASK,res2)
 

  prcMASK = where(lsmTRMM.eq.0, TRMM_reorder(nt,:,:), TRMM_reorder@_FillValue) ; Ocean only
  plot(2) = gsn_csm_contour_map(wks,prcMASK,res2)
  
  panres = True                                ; create a panel resource list
  panres@gsnMaximize = True                    ; maximize the size of the paneled plots
  panres@gsnPanelYWhiteSpacePercent = 1.0      ; leave space between the 2 plots
  panres@gsnPanelLabelBar = True               ; add common colorbar
 ;panres@lbLabelFontHeightF  = 0.007             ; make labels smaller
  gsn_panel(wks,plot,(/3,1/),panres)
end if