; ***********************************************

 ; 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