npp/03.0.05deg_to_T31.ncl
changeset 0 0c6405ab2ff4
equal deleted inserted replaced
-1:000000000000 0:a4abee0712e9
       
     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 ; output data
       
    11 ;************************************************
       
    12   diro  = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"
       
    13   filo  = "Npp_T31_mean.nc"
       
    14   c = addfile(diro+filo,"c")
       
    15 
       
    16 ;************************************************
       
    17 ; read in observed data
       
    18 ;************************************************
       
    19   diri  = "/fis/cgd/cseg/people/jeff/clamp_data/npp/ob/"
       
    20   fili  = "Npp_0.05deg_mean.nc"
       
    21   g     = addfile (diri+fili,"r")
       
    22   bi    = g->NPP   
       
    23   xi    = g->lon 
       
    24   yi    = g->lat
       
    25 
       
    26 ;************************************************
       
    27 ; change into 0-360E, 90S-90N
       
    28 ; Observed NPP*scale_factor
       
    29 ;************************************************
       
    30 
       
    31  scale_factor = 0.1
       
    32  
       
    33  yi    = (/ yi(::-1) /)
       
    34  bi    =  (/ bi(::-1,:) /)
       
    35  b2    = bi
       
    36  x2    = xi   
       
    37  
       
    38  nx = dimsizes(xi)
       
    39  do i= 0,nx-1
       
    40     if (i .lt. 3600) then
       
    41        p = i + 3600
       
    42        xi(p) = x2(i) + 360.      
       
    43     else
       
    44        p = i - 3600
       
    45        xi(p) = x2(i)
       
    46     end if
       
    47     bi(:,p)= b2(:,i) * scale_factor
       
    48  end do
       
    49 
       
    50  bi&lat =  yi
       
    51  bi&lon =  xi
       
    52 
       
    53 ;print (xi)
       
    54 ;print (yi)
       
    55 ;exit
       
    56 
       
    57 ;************************************************
       
    58 ; read in model data
       
    59 ;************************************************
       
    60  diri2  = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
       
    61  fili2  = "b30.061n_1995-2004_ANN_climo_lnd.nc"
       
    62  f      = addfile (diri2+fili2,"r")
       
    63 
       
    64  lon    = f->lon     
       
    65  lat    = f->lat
       
    66  nlon   = dimsizes(lon)
       
    67  nlat   = dimsizes(lat)      
       
    68 
       
    69 ; print (xi)
       
    70 ; print (yi)
       
    71 ; print (xo)
       
    72 ; print (yo)
       
    73 
       
    74 ;(A)
       
    75 ;bo = linint2_Wrap(xi,yi,bi,True,xo,yo,0)
       
    76 
       
    77 ;(B)
       
    78 
       
    79   rad   = 4.*atan(1.)/180.
       
    80   clat  = lat
       
    81   clat  = cos(lat*rad)
       
    82   clat@long_name = "cos(latitude)"
       
    83   delete(clat@units)
       
    84   printVarSummary(clat)
       
    85 
       
    86  bo = new((/nlat,nlon/),float)
       
    87     do j=0,nlat-1
       
    88        if (j.eq.0 .or. j.eq.nlat-1) then
       
    89           if (j.eq.0) then
       
    90              LATS = -90.          
       
    91              LATN = lat(j)+0.5*(lat(j+1)-lat(j))
       
    92           end if
       
    93           if (j.eq.nlat-1) then
       
    94              LATS = lat(j)-0.5*(lat(j)-lat(j-1))
       
    95              LATN = 90.                  
       
    96           end if
       
    97        else
       
    98           LATS = lat(j)-0.5*(lat(j)-lat(j-1))
       
    99           LATN = lat(j)+0.5*(lat(j+1)-lat(j))
       
   100        end if
       
   101  
       
   102 ;      CLAT = clat({LATS:LATN})      ; do once for *slight* efficiency
       
   103 ;      TEMP = bi({LATS:LATN},:)      ; 2D [lat,lon]
       
   104  
       
   105       do i=0,nlon-1
       
   106        if (i.eq.0 .or. i.eq.nlon-1) then
       
   107           if (i.eq.0) then
       
   108              LONL = 0.          
       
   109              LONR = lon(i)+0.5*(lon(i+1)-lon(i))
       
   110           end if
       
   111           if (i.eq.nlon-1) then
       
   112              LONL = lon(i)-0.5*(lon(i)-lon(i-1))
       
   113              LONR = 360.                 
       
   114           end if
       
   115        else
       
   116           LONL = lon(i)-0.5*(lon(i)-lon(i-1))
       
   117           LONR = lon(i)+0.5*(lon(i+1)-lon(i))
       
   118        end if
       
   119 
       
   120 ;print (LATS)
       
   121 ;print (LATN)
       
   122 ;print (LONL)
       
   123 ;print (LONR)
       
   124          bo(j,i) = avg(bi({LATS:LATN},{LONL:LONR}))  
       
   125 ;        bo(j,i) = wgt_areaave(TEMP(:,{LONL:LONR}), CLAT, 1.0, 0)
       
   126       end do
       
   127  
       
   128 ;     delete(CLAT)
       
   129 ;     delete(TEMP) 
       
   130     end do
       
   131 
       
   132   bo!0   = "lat"
       
   133   bo!1   = "lon"
       
   134   bo&lat =  lat
       
   135   bo&lon =  lon
       
   136   bo@units      = bi@units
       
   137   bo@long_name  = bi@long_name
       
   138 ; bo@_FillValue = bi@_FillValue
       
   139   bo@_FillValue = 1.e+36
       
   140 
       
   141   c->NPP  = bo
       
   142 end