load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"

begin

 lat=new((/360/),float)
 lon=new((/720/),float)
 val=new((/1/),float)
 val=90.
 do ilat=0,359
   lat(ilat)=val
   val=val-0.5
 end do
 val=-180.
 do ilon=0,719
   lon(ilon)=val
   val=val+0.5
 end do
 lat@units="degrees_north"
 lon@units="degrees_east"

 data = asciiread("/ptmp/notaro/potveg/potential_veg_hd.asc",(/360,720/),"float")
 data!0="lat"
 data!1="lon"
 data&lat=lat
 data&lon=lon
 data@_FillValue=0

  biome=new((/360,720/),float)
  do ilat=0,359
    do ilon=0,719
      if (.not.ismissing(data(ilat,ilon))) then
        if ((data(ilat,ilon).eq.1).or.(data(ilat,ilon).eq.2)) then
          biome(ilat,ilon)=1. ; tropical forest
        end if
        if ((data(ilat,ilon).ge.3).and.(data(ilat,ilon).le.5)) then
          biome(ilat,ilon)=2. ; temperate forest
        end if
        if ((data(ilat,ilon).eq.6).or.(data(ilat,ilon).eq.7)) then
          biome(ilat,ilon)=3. ; boreal forest
        end if
        if (data(ilat,ilon).eq.9) then
          biome(ilat,ilon)=4. ; savanna/dry woodland
        end if
        if ((data(ilat,ilon).ge.10).and.(data(ilat,ilon).le.12)) then
          biome(ilat,ilon)=5. ; dry grassland/shrubland
        end if
        if (data(ilat,ilon).eq.14) then
          biome(ilat,ilon)=6. ; desert
        end if
        if (data(ilat,ilon).eq.13) then
          biome(ilat,ilon)=7. ; tundra
        end if
;        if (data(ilat,ilon).eq.10) then
;          biome(ilat,ilon)=1e+35 ; polar desert
;        end if

        if (data(ilat,ilon).eq.8) then
          if ((lat(ilat).lt.25.) .and. (lat(ilat).gt.-25.)) then
            biome(ilat,ilon)=1.
          end if
          if ((fabs(lat(ilat)).ge.25.) .and. (fabs(lat(ilat)).le.50.)) then
            biome(ilat,ilon)=2.
          end if
          if (fabs(lat(ilat)).gt.50.) then
            biome(ilat,ilon)=3.
          end if
        end if

      else
        biome(ilat,ilon)=1e+35
      end if
    end do
  end do

  biome!0="lat"
  biome!1="lon"
  biome&lat=lat
  biome&lon=lon
  biome@_FillValue=1e+35

  wks = gsn_open_wks("ps","biomes_potveg_plot_recolor")          ; open a ps file
  res                  = True                ; plot mods desired
  res@cnFillOn         = True               ; color Fill
  res@cnRasterModeOn   = True               ; Raster Mode
  res@cnLinesOn        =  False             ; Turn off contour lines
  res@tmXTOn=False
  res@tmYROn=False
  res@cnLevelSelectionMode = "ManualLevels" ; set manual contour levels
  res@cnMinLevelValF   = 1.0               ; set min contour level
  res@cnMaxLevelValF   =  6.0               ; set max contour level
  res@cnLevelSpacingF  =  1.0               ; set contour spacing
  res@gsnCenterString  = "Potential Vegetation Biome Map"   ; add center string
  res@mpMinLatF = -60.
  res@mpMaxLatF = 75.
;;  colors = \
;;   (/"white","black","pale green","orange","dark green","blue","red","yellow","grey"/)
  colors = \
   (/"white","black","pale green","forest green","royal blue","red","orange","yellow","grey"/)
   gsn_define_colormap(wks,colors)
   res@lbLabelStrings =  (/"","","","","","","",""/)
  txres = True   ; Set up variable for TextItem resources.
  txres@txFontHeightF = 0.019                ; Set the font height.
  gsn_text_ndc(wks,"1",.23,.25,txres)
  gsn_text_ndc(wks,"2",.32,.25,txres)
  gsn_text_ndc(wks,"3",.41,.25,txres)
  gsn_text_ndc(wks,"4",.5,.25,txres)
  gsn_text_ndc(wks,"5",.59,.25,txres)
  gsn_text_ndc(wks,"6",.68,.25,txres)
  gsn_text_ndc(wks,"7",.77,.25,txres)
  txres@txFontHeightF = 0.01                ; Set the font height.
  txres@txJust = "CenterLeft"
  gsn_text_ndc(wks,"Cat 1 - Tropical forest",.08,.20,txres)
  gsn_text_ndc(wks,"Cat 2 - Temperate forest",.08,.17,txres)
  gsn_text_ndc(wks,"Cat 3 - Boreal forest",.08,.14,txres)
  gsn_text_ndc(wks,"Cat 4 - Savanna / Dry woodland",.08,.11,txres)
  gsn_text_ndc(wks,"Cat 5 - Dry grassland / Shrubland",.56,.2,txres)
  gsn_text_ndc(wks,"Cat 6 - Desert",.56,.17,txres)
  gsn_text_ndc(wks,"Cat 7 - Tundra",.56,.14,txres)

  plot = gsn_csm_contour_map_ce(wks,biome, res) ; create plot

  end