lai/03.lai_ob_to_T31.ncl.x
changeset 1 4be95183fbcd
equal deleted inserted replaced
-1:000000000000 0:d50cd7ef2d00
       
     1 ; ***********************************************
       
     2 ; interpolate into model grids (T31)
       
     3 ; ***********************************************
       
     4 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
       
     5 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
       
     6 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
       
     7 ;************************************************
       
     8 begin
       
     9 
       
    10 ;************************************************
       
    11 ; output data
       
    12 ;************************************************
       
    13   diro  = "/fis/cgd/cseg/people/jeff/clamp_data/lai/ob/"
       
    14   filo  = "LAI_2000-2005_ensemble_T31.nc"
       
    15   c = addfile(diro+filo,"c")
       
    16   filedimdef(c,"time",-1,True)
       
    17 
       
    18 ;************************************************
       
    19 ; read in observed data
       
    20 ;************************************************
       
    21   diri  = "/fis/cgd/cseg/people/jeff/clamp_data/lai/ob/"
       
    22   fili  = "LAI_2000-2005_ensemble_0.05deg.nc"
       
    23   g     = addfile (diri+fili,"r")
       
    24   bi    = g->LAI   
       
    25   xi    = g->lon 
       
    26   yi    = g->lat
       
    27 
       
    28 ;************************************************
       
    29 ; change into 0-360E, 90S-90N
       
    30 ;************************************************
       
    31  
       
    32  yi    = (/ yi(::-1) /)
       
    33  bi    = (/ bi(:,::-1,:) /)
       
    34  printVarSummary(bi)
       
    35 
       
    36  b2    = bi
       
    37  x2    = xi   
       
    38  
       
    39  nx = dimsizes(xi)
       
    40  do i= 0,nx-1
       
    41     if (i .lt. 3600) then
       
    42        p = i + 3600
       
    43        xi(p) = x2(i) + 360.      
       
    44     else
       
    45        p = i - 3600
       
    46        xi(p) = x2(i)
       
    47     end if
       
    48     bi(:,:,p)= b2(:,:,i) 
       
    49  end do
       
    50 
       
    51  bi&lat =  yi
       
    52  bi&lon =  xi
       
    53 
       
    54 ;print (xi)
       
    55 ;print (yi)
       
    56 ;************************************************
       
    57 ; read in model data
       
    58 ;************************************************
       
    59  diri2  = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
       
    60  fili2  = "b30.061n_1995-2004_MONS_climo_lnd.nc"
       
    61  f     = addfile (diri2+fili2,"r")
       
    62 
       
    63  lon    = f->lon     
       
    64  lat    = f->lat
       
    65  nlon   = dimsizes(lon)
       
    66  nlat   = dimsizes(lat)      
       
    67 
       
    68 ; print (xi)
       
    69 ; print (yi)
       
    70 ; print (xo)
       
    71 ; print (yo)
       
    72 
       
    73  bo = new((/12,nlat,nlon/),float)
       
    74 
       
    75  do m = 0,11
       
    76     do j=0,nlat-1
       
    77        if (j.eq.0 .or. j.eq.nlat-1) then
       
    78           if (j.eq.0) then
       
    79              LATS = -90.          
       
    80              LATN = lat(j)+0.5*(lat(j+1)-lat(j))
       
    81           end if
       
    82           if (j.eq.nlat-1) then
       
    83              LATS = lat(j)-0.5*(lat(j)-lat(j-1))
       
    84              LATN = 90.                  
       
    85           end if
       
    86        else
       
    87           LATS = lat(j)-0.5*(lat(j)-lat(j-1))
       
    88           LATN = lat(j)+0.5*(lat(j+1)-lat(j))
       
    89        end if
       
    90  
       
    91       do i=0,nlon-1
       
    92        if (i.eq.0 .or. i.eq.nlon-1) then
       
    93           if (i.eq.0) then
       
    94              LONL = 0.          
       
    95              LONR = lon(i)+0.5*(lon(i+1)-lon(i))
       
    96           end if
       
    97           if (i.eq.nlon-1) then
       
    98              LONL = lon(i)-0.5*(lon(i)-lon(i-1))
       
    99              LONR = 360.                 
       
   100           end if
       
   101        else
       
   102           LONL = lon(i)-0.5*(lon(i)-lon(i-1))
       
   103           LONR = lon(i)+0.5*(lon(i+1)-lon(i))
       
   104        end if
       
   105 
       
   106 ;print (LATS)
       
   107 ;print (LATN)
       
   108 ;print (LONL)
       
   109 ;print (LONR)
       
   110          bo(m,j,i) = avg(bi(m,{LATS:LATN},{LONL:LONR}))  
       
   111       end do 
       
   112     end do
       
   113   end do
       
   114 
       
   115   bo!0   = "time"
       
   116   bo!1   = "lat"
       
   117   bo!2   = "lon"
       
   118   bo&time= bi&time
       
   119   bo&lat = lat
       
   120   bo&lon = lon
       
   121   bo@units      = bi@units
       
   122   bo@long_name  = bi@long_name
       
   123   bo@_FillValue = bi@_FillValue
       
   124 ; bo@units      = "none"
       
   125 ; bo@long_name  = "Leaf Area Index"
       
   126 
       
   127   c->LAI  = bo
       
   128 end