1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/lai/01.read.land_class.split.ncl	Mon Jan 26 22:08:20 2009 -0500
     1.3 @@ -0,0 +1,67 @@
     1.4 +;************************************************
     1.5 +;    Read ascii, Write nc
     1.6 +;    output: lat: N->S     lon: -180W->180E 
     1.7 +;    class 8 split into class  8 = woody savanna (S. Hem.)
     1.8 +;                             18 = woody savanna (N. Hem.)                   
     1.9 +;    class 8 split into class  9 =       savanna (S. Hem.)
    1.10 +;                             19 =       savanna (N. Hem.)     
    1.11 +;************************************************
    1.12 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
    1.13 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
    1.14 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"  
    1.15 +;************************************************
    1.16 +begin
    1.17 +;---------------------------------------------------
    1.18 +; final data
    1.19 +  diro  = "/fis/cgd/cseg/people/jeff/clamp_data/lai/"
    1.20 +  filo  = "land_class_0.25deg_new.nc"
    1.21 +  c = addfile(diro+filo,"c")
    1.22 +; filedimdef(c,"time",-1,True) 
    1.23 +       
    1.24 +  nlat  = 180*4 
    1.25 +  nlon  = 360*4
    1.26 +  
    1.27 +  t        = new((/nlat,nlon/),integer)
    1.28 +  lon      = new((/nlon/),float)
    1.29 +  lat      = new((/nlat/),float)
    1.30 +;---------------------------------------------------
    1.31 +; input data
    1.32 +  diri  = "/fis/cgd/cseg/people/jeff/clamp_data/lai/modis_landcover_qdeg/"
    1.33 +  fili = "modis_landcover_class_qd.asc"                    
    1.34 +  b = diri+fili
    1.35 +  
    1.36 +  lat  = latGlobeFo(nlat, "lat", "latitude", "degrees_north")
    1.37 +  lat  = (/ lat(::-1) /)                     ; make N->S
    1.38 +  lon  = lonGlobeFo(nlon, "lon", "longitude", "degrees_east")
    1.39 +  lon  = (/ lon - 180. /) ; subtract 180 from all values 
    1.40 +  lon&lon = lon           ; update coordinates
    1.41 +  lat&lat = lat           ; update coordinates
    1.42 +
    1.43 +;=============================
    1.44 +;  name dimensions of t and assign coordinate variables
    1.45 +;============================  
    1.46 +  t!0    = "lat"
    1.47 +  t!1    = "lon"
    1.48 +  t&lat  = lat
    1.49 +  t&lon  = lon
    1.50 +  t@long_name = "landcover class"
    1.51 +  t@units     = ""
    1.52 +  t@_FillValue= -9999         
    1.53 +  t@missing_value= -9999         
    1.54 +
    1.55 +  t = asciiread(b,(/nlat,nlon/),"integer")
    1.56 +
    1.57 +  nlat = dimsizes(lat)
    1.58 +  ny2  = nlat/2
    1.59 +
    1.60 +  do j = 0,ny2-1
    1.61 +     t(j,:) = where((t(j,:) .eq. 8),18,t(j,:))
    1.62 +     t(j,:) = where((t(j,:) .eq. 9),19,t(j,:)) 
    1.63 +  end do
    1.64 + 
    1.65 +  c->LAND_CLASS  = t
    1.66 +  c->lat  = lat
    1.67 +  c->lon  = lon
    1.68 +
    1.69 +  print (min(t) + "/" + max(t))
    1.70 +end