lai/03.land_class_to_T42.ncl
changeset 1 4be95183fbcd
equal deleted inserted replaced
-1:000000000000 0:4b1b19104037
       
     1 ; ***********************************************
       
     2 ; interpolate into model grids (T42)
       
     3 ; output: lat: S->N    lon: 0->360
       
     4 ; input : lat: N->S    lon: -180->180 
       
     5 ; ***********************************************
       
     6 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
       
     7 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
       
     8 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
       
     9 ;************************************************
       
    10 begin
       
    11 
       
    12   nclass = 20
       
    13 
       
    14 ;************************************************
       
    15 ; output data
       
    16 ;************************************************
       
    17   diro  = "/fis/cgd/cseg/people/jeff/clamp_data/lai/"
       
    18 ; filo  = "land_class_T42.nc"
       
    19   filo  = "land_class_T42_new.nc"
       
    20   c = addfile(diro+filo,"c")
       
    21 
       
    22 ;************************************************
       
    23 ; read in observed data
       
    24 ;************************************************
       
    25   diri  = "/fis/cgd/cseg/people/jeff/clamp_data/lai/"
       
    26   fili  = "land_class_0.25deg.nc"
       
    27   fili  = "land_class_0.25deg_new.nc"
       
    28   g     = addfile (diri+fili,"r")
       
    29   bi    = g->LAND_CLASS   
       
    30   xi    = g->lon 
       
    31   yi    = g->lat
       
    32 
       
    33 ;************************************************
       
    34 ; change into 0-360E, 90S-90N
       
    35 ;************************************************
       
    36  
       
    37  yi    = (/ yi(::-1) /)
       
    38  bi    = (/ bi(::-1,:) /)
       
    39  printVarSummary(bi)
       
    40 
       
    41  b2    = bi
       
    42  x2    = xi   
       
    43  
       
    44  nx  = dimsizes(xi)
       
    45  nx2 = nx/2
       
    46 
       
    47  do i= 0,nx-1
       
    48     if (i .lt. nx2) then
       
    49        p = i + nx2
       
    50        xi(p) = x2(i) + 360.      
       
    51     else
       
    52        p = i - nx2
       
    53        xi(p) = x2(i)
       
    54     end if
       
    55     bi(:,p)= b2(:,i) 
       
    56  end do
       
    57 
       
    58  bi&lat =  yi
       
    59  bi&lon =  xi
       
    60 
       
    61  print (xi)
       
    62  print (yi)
       
    63 ;exit
       
    64 
       
    65 ;************************************************
       
    66 ; read in model data
       
    67 ;************************************************
       
    68  diri2  = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
       
    69  fili2  = "i01.03cn_1545-1569_ANN_climo.nc"
       
    70  f     = addfile (diri2+fili2,"r")
       
    71 
       
    72  lon    = f->lon     
       
    73  lat    = f->lat
       
    74  nlon   = dimsizes(lon)
       
    75  nlat   = dimsizes(lat)      
       
    76 
       
    77 ; print (xi)
       
    78 ; print (yi)
       
    79 ; print (xo)
       
    80 ; print (yo)
       
    81 
       
    82  bo = new((/nlat,nlon/),integer)
       
    83  count = new((/nclass/),integer)
       
    84 
       
    85   do j=0,nlat-1
       
    86      if (j.eq.0 .or. j.eq.nlat-1) then
       
    87         if (j.eq.0) then
       
    88            LATS = -90.          
       
    89            LATN = lat(j)+0.5*(lat(j+1)-lat(j))
       
    90         end if
       
    91         if (j.eq.nlat-1) then
       
    92            LATS = lat(j)-0.5*(lat(j)-lat(j-1))
       
    93            LATN = 90.                  
       
    94         end if
       
    95      else
       
    96         LATS = lat(j)-0.5*(lat(j)-lat(j-1))
       
    97         LATN = lat(j)+0.5*(lat(j+1)-lat(j))
       
    98      end if
       
    99  
       
   100      do i=0,nlon-1
       
   101         if (i.eq.0 .or. i.eq.nlon-1) then
       
   102            if (i.eq.0) then
       
   103               LONL = 0.          
       
   104               LONR = lon(i)+0.5*(lon(i+1)-lon(i))
       
   105            end if
       
   106            if (i.eq.nlon-1) then
       
   107               LONL = lon(i)-0.5*(lon(i)-lon(i-1))
       
   108               LONR = 360.                 
       
   109            end if
       
   110         else
       
   111            LONL = lon(i)-0.5*(lon(i)-lon(i-1))
       
   112            LONR = lon(i)+0.5*(lon(i+1)-lon(i))
       
   113         end if
       
   114 
       
   115 ;print (LATS)
       
   116 ;print (LATN)
       
   117 ;print (LONL)
       
   118 ;print (LONR)
       
   119 
       
   120         count = 0
       
   121 
       
   122         do n = 0,nclass-1
       
   123            count(n) = num (bi({LATS:LATN},{LONL:LONR}).eq.n)
       
   124         end do
       
   125  
       
   126         bo(j,i) = maxind(count)  
       
   127      end do     
       
   128   end do
       
   129 
       
   130   bo!0   = "lat"
       
   131   bo!1   = "lon"
       
   132   bo&lat = lat
       
   133   bo&lon = lon
       
   134   bo@units      = bi@units
       
   135   bo@long_name  = bi@long_name
       
   136   bo@_FillValue = bi@_FillValue
       
   137 
       
   138   print (min(bo) + "/" + max(bo))
       
   139 
       
   140   c->LAND_CLASS  = bo
       
   141 end