forrest@0: ;************************************************
forrest@0: ;    Read ascii, Write nc
forrest@0: ;    output: lat: N->S     lon: -180W->180E 
forrest@0: ;    class 8 split into class  8 = woody savanna (S. Hem.)
forrest@0: ;                             18 = woody savanna (N. Hem.)                   
forrest@0: ;    class 8 split into class  9 =       savanna (S. Hem.)
forrest@0: ;                             19 =       savanna (N. Hem.)     
forrest@0: ;************************************************
forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
forrest@0: load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"  
forrest@0: ;************************************************
forrest@0: begin
forrest@0: ;---------------------------------------------------
forrest@0: ; final data
forrest@0:   diro  = "/fis/cgd/cseg/people/jeff/clamp_data/lai/"
forrest@0:   filo  = "land_class_0.25deg_new.nc"
forrest@0:   c = addfile(diro+filo,"c")
forrest@0: ; filedimdef(c,"time",-1,True) 
forrest@0:        
forrest@0:   nlat  = 180*4 
forrest@0:   nlon  = 360*4
forrest@0:   
forrest@0:   t        = new((/nlat,nlon/),integer)
forrest@0:   lon      = new((/nlon/),float)
forrest@0:   lat      = new((/nlat/),float)
forrest@0: ;---------------------------------------------------
forrest@0: ; input data
forrest@0:   diri  = "/fis/cgd/cseg/people/jeff/clamp_data/lai/modis_landcover_qdeg/"
forrest@0:   fili = "modis_landcover_class_qd.asc"                    
forrest@0:   b = diri+fili
forrest@0:   
forrest@0:   lat  = latGlobeFo(nlat, "lat", "latitude", "degrees_north")
forrest@0:   lat  = (/ lat(::-1) /)                     ; make N->S
forrest@0:   lon  = lonGlobeFo(nlon, "lon", "longitude", "degrees_east")
forrest@0:   lon  = (/ lon - 180. /) ; subtract 180 from all values 
forrest@0:   lon&lon = lon           ; update coordinates
forrest@0:   lat&lat = lat           ; update coordinates
forrest@0: 
forrest@0: ;=============================
forrest@0: ;  name dimensions of t and assign coordinate variables
forrest@0: ;============================  
forrest@0:   t!0    = "lat"
forrest@0:   t!1    = "lon"
forrest@0:   t&lat  = lat
forrest@0:   t&lon  = lon
forrest@0:   t@long_name = "landcover class"
forrest@0:   t@units     = ""
forrest@0:   t@_FillValue= -9999         
forrest@0:   t@missing_value= -9999         
forrest@0: 
forrest@0:   t = asciiread(b,(/nlat,nlon/),"integer")
forrest@0: 
forrest@0:   nlat = dimsizes(lat)
forrest@0:   ny2  = nlat/2
forrest@0: 
forrest@0:   do j = 0,ny2-1
forrest@0:      t(j,:) = where((t(j,:) .eq. 8),18,t(j,:))
forrest@0:      t(j,:) = where((t(j,:) .eq. 9),19,t(j,:)) 
forrest@0:   end do
forrest@0:  
forrest@0:   c->LAND_CLASS  = t
forrest@0:   c->lat  = lat
forrest@0:   c->lon  = lon
forrest@0: 
forrest@0:   print (min(t) + "/" + max(t))
forrest@0: end