biomass/03.ncl
changeset 1 4be95183fbcd
equal deleted inserted replaced
-1:000000000000 0:64580d5834e7
       
     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 ;************************************************
       
    12 ; read in observed data
       
    13 ;************************************************
       
    14   diri  = "/fis/cgd/cseg/people/jeff/clamp_data/biomass/ob/"
       
    15   fili  = "amazon_biomass.nc"
       
    16   g     = addfile (diri+fili,"r")
       
    17   bi    = g->BIOMASS   
       
    18   xi    = g->lon 
       
    19   yi    = g->lat
       
    20 
       
    21 
       
    22 ;************************************************
       
    23 ; read in model data
       
    24 ;************************************************
       
    25  diri2  = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
       
    26  fili2  = "b30.061n_1995-2004_ANN_climo_lnd.nc"
       
    27  f      = addfile (diri2+fili2,"r")
       
    28 
       
    29  lon    = f->lon     
       
    30  lat    = f->lat
       
    31  nlon   = dimsizes(lon)
       
    32  nlat   = dimsizes(lat)      
       
    33 
       
    34 ; print (xi)
       
    35 ; print (yi)
       
    36 
       
    37  ind_lonL = min(ind(lon .ge. xi(0)))
       
    38  ind_lonR = max(ind(lon .le. xi(dimsizes(xi)-1)))
       
    39 
       
    40 print (ind_lonL)
       
    41 print (xi(0))
       
    42 print (lon(ind_lonL))
       
    43 
       
    44 print (ind_lonR)
       
    45 print (xi(dimsizes(xi)-1))
       
    46 print (lon(ind_lonR))
       
    47 
       
    48  ind_latS = min(ind(lat .ge. yi(0)))
       
    49  ind_latN = max(ind(lat .le. yi(dimsizes(yi)-1)))
       
    50 
       
    51 print (ind_latS)
       
    52 print (yi(0))
       
    53 print (lat(ind_latS))
       
    54 
       
    55 print (ind_latN)
       
    56 print (yi(dimsizes(yi)-1))
       
    57 print (lat(ind_latN))
       
    58 exit
       
    59  bo = new((/nlat,nlon/),float)
       
    60     do j=0,nlat-1
       
    61        if (j.eq.0 .or. j.eq.nlat-1) then
       
    62           if (j.eq.0) then
       
    63              LATS = -90.          
       
    64              LATN = lat(j)+0.5*(lat(j+1)-lat(j))
       
    65           end if
       
    66           if (j.eq.nlat-1) then
       
    67              LATS = lat(j)-0.5*(lat(j)-lat(j-1))
       
    68              LATN = 90.                  
       
    69           end if
       
    70        else
       
    71           LATS = lat(j)-0.5*(lat(j)-lat(j-1))
       
    72           LATN = lat(j)+0.5*(lat(j+1)-lat(j))
       
    73        end if
       
    74  
       
    75 ;      CLAT = clat({LATS:LATN})      ; do once for *slight* efficiency
       
    76 ;      TEMP = bi({LATS:LATN},:)      ; 2D [lat,lon]
       
    77  
       
    78       do i=0,nlon-1
       
    79        if (i.eq.0 .or. i.eq.nlon-1) then
       
    80           if (i.eq.0) then
       
    81              LONL = 0.          
       
    82              LONR = lon(i)+0.5*(lon(i+1)-lon(i))
       
    83           end if
       
    84           if (i.eq.nlon-1) then
       
    85              LONL = lon(i)-0.5*(lon(i)-lon(i-1))
       
    86              LONR = 360.                 
       
    87           end if
       
    88        else
       
    89           LONL = lon(i)-0.5*(lon(i)-lon(i-1))
       
    90           LONR = lon(i)+0.5*(lon(i+1)-lon(i))
       
    91        end if
       
    92 
       
    93 ;print (LATS)
       
    94 ;print (LATN)
       
    95 ;print (LONL)
       
    96 ;print (LONR)
       
    97          bo(j,i) = avg(bi({LATS:LATN},{LONL:LONR}))  
       
    98 ;        bo(j,i) = wgt_areaave(TEMP(:,{LONL:LONR}), CLAT, 1.0, 0)
       
    99       end do
       
   100  
       
   101 ;     delete(CLAT)
       
   102 ;     delete(TEMP) 
       
   103     end do
       
   104 
       
   105   bo!0   = "lat"
       
   106   bo!1   = "lon"
       
   107   bo&lat =  lat
       
   108   bo&lon =  lon
       
   109   bo@units      = bi@units
       
   110   bo@long_name  = bi@long_name
       
   111 ; bo@_FillValue = bi@_FillValue
       
   112   bo@_FillValue = 1.e+36
       
   113 
       
   114   c->NPP  = bo
       
   115 end