diff -r 000000000000 -r 0c6405ab2ff4 lai/03.land_class_to_T31.ncl
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/lai/03.land_class_to_T31.ncl	Mon Jan 26 22:08:20 2009 -0500
@@ -0,0 +1,139 @@
+; ***********************************************
+; interpolate into model grids (T42)
+; output: lat: S->N    lon: 0->360
+; input : lat: N->S    lon: -180->180 
+; ***********************************************
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
+;************************************************
+begin
+
+  nclass = 20
+
+;************************************************
+; output data
+;************************************************
+  diro  = "/fis/cgd/cseg/people/jeff/clamp_data/lai/ob/"
+  filo  = "land_class_T31.nc"
+  c = addfile(diro+filo,"c")
+
+;************************************************
+; read in observed data
+;************************************************
+  diri  = "/fis/cgd/cseg/people/jeff/clamp_data/lai/ob/"
+  fili  = "land_class_0.25deg_new.nc"
+  g     = addfile (diri+fili,"r")
+  bi    = g->LAND_CLASS   
+  xi    = g->lon 
+  yi    = g->lat
+
+;************************************************
+; change into 0-360E, 90S-90N
+;************************************************
+ 
+ yi    = (/ yi(::-1) /)
+ bi    = (/ bi(::-1,:) /)
+ printVarSummary(bi)
+
+ b2    = bi
+ x2    = xi   
+ 
+ nx  = dimsizes(xi)
+ nx2 = nx/2
+
+ do i= 0,nx-1
+    if (i .lt. nx2) then
+       p = i + nx2
+       xi(p) = x2(i) + 360.      
+    else
+       p = i - nx2
+       xi(p) = x2(i)
+    end if
+    bi(:,p)= b2(:,i) 
+ end do
+
+ bi&lat =  yi
+ bi&lon =  xi
+
+ print (xi)
+ print (yi)
+;exit
+
+;************************************************
+; read in model data
+;************************************************
+ diri2  = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
+ fili2  = "b30.061n_1995-2004_ANN_climo_lnd.nc"
+ f     = addfile (diri2+fili2,"r")
+
+ lon    = f->lon     
+ lat    = f->lat
+ nlon   = dimsizes(lon)
+ nlat   = dimsizes(lat)      
+
+; print (xi)
+; print (yi)
+; print (xo)
+; print (yo)
+
+ bo = new((/nlat,nlon/),integer)
+ count = new((/nclass/),integer)
+
+  do j=0,nlat-1
+     if (j.eq.0 .or. j.eq.nlat-1) then
+        if (j.eq.0) then
+           LATS = -90.          
+           LATN = lat(j)+0.5*(lat(j+1)-lat(j))
+        end if
+        if (j.eq.nlat-1) then
+           LATS = lat(j)-0.5*(lat(j)-lat(j-1))
+           LATN = 90.                  
+        end if
+     else
+        LATS = lat(j)-0.5*(lat(j)-lat(j-1))
+        LATN = lat(j)+0.5*(lat(j+1)-lat(j))
+     end if
+ 
+     do i=0,nlon-1
+        if (i.eq.0 .or. i.eq.nlon-1) then
+           if (i.eq.0) then
+              LONL = 0.          
+              LONR = lon(i)+0.5*(lon(i+1)-lon(i))
+           end if
+           if (i.eq.nlon-1) then
+              LONL = lon(i)-0.5*(lon(i)-lon(i-1))
+              LONR = 360.                 
+           end if
+        else
+           LONL = lon(i)-0.5*(lon(i)-lon(i-1))
+           LONR = lon(i)+0.5*(lon(i+1)-lon(i))
+        end if
+
+;print (LATS)
+;print (LATN)
+;print (LONL)
+;print (LONR)
+
+        count = 0
+
+        do n = 0,nclass-1
+           count(n) = num (bi({LATS:LATN},{LONL:LONR}).eq.n)
+        end do
+ 
+        bo(j,i) = maxind(count)  
+     end do     
+  end do
+
+  bo!0   = "lat"
+  bo!1   = "lon"
+  bo&lat = lat
+  bo&lon = lon
+  bo@units      = bi@units
+  bo@long_name  = bi@long_name
+  bo@_FillValue = bi@_FillValue
+
+  print (min(bo) + "/" + max(bo))
+
+  c->LAND_CLASS  = bo
+end
\ No newline at end of file