lai/03.lai_ob_to_T42.ncl
changeset 1 4be95183fbcd
equal deleted inserted replaced
-1:000000000000 0:94f6ba34ef5d
       
     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   year = 2003
       
    11 
       
    12 ;************************************************
       
    13 ; output data
       
    14 ;************************************************
       
    15   diro  = "/fis/cgd/cseg/people/jeff/clamp_data/"
       
    16 ; filo  = "LAI_" + year + "_monthly_T42.nc"
       
    17 ; filo  = "LAI_2000-2005_mean_T42.nc"
       
    18   filo  = "LAI_2000-2005_ensemble_T42.nc"
       
    19   c = addfile(diro+filo,"c")
       
    20   filedimdef(c,"time",-1,True)
       
    21 
       
    22 ;************************************************
       
    23 ; read in observed data
       
    24 ;************************************************
       
    25   diri  = "/fis/cgd/cseg/people/jeff/clamp_data/"
       
    26 ; fili  = "LAI_" + year + "_monthly.nc"
       
    27 ; fili  = "LAI_2000-2005_mean_0.05deg.nc"
       
    28   fili  = "LAI_2000-2005_ensemble_0.05deg.nc"
       
    29   g     = addfile (diri+fili,"r")
       
    30   bi    = g->LAI   
       
    31   xi    = g->lon 
       
    32   yi    = g->lat
       
    33 
       
    34 ;************************************************
       
    35 ; change into 0-360E, 90S-90N
       
    36 ; Observed NPP*scale_factor
       
    37 ;************************************************
       
    38  
       
    39  yi    = (/ yi(::-1) /)
       
    40  bi    = (/ bi(:,::-1,:) /)
       
    41  printVarSummary(bi)
       
    42 
       
    43  b2    = bi
       
    44  x2    = xi   
       
    45  
       
    46  nx = dimsizes(xi)
       
    47  do i= 0,nx-1
       
    48     if (i .lt. 3600) then
       
    49        p = i + 3600
       
    50        xi(p) = x2(i) + 360.      
       
    51     else
       
    52        p = i - 3600
       
    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 ;fili2  = "i01.04casa_1605-1629_ANN_climo.nc"
       
    69  fili2  = "i01.03cn_1545-1569_ANN_climo.nc"
       
    70  f     = addfile (diri+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 ;(A)
       
    83 ;bo = linint2_Wrap(xi,yi,bi,True,xo,yo,0)
       
    84 
       
    85 ;(B)
       
    86 
       
    87 ; rad   = 4.*atan(1.)/180.
       
    88 ; clat  = lat
       
    89 ; clat  = cos(lat*rad)
       
    90 ; clat@long_name = "cos(latitude)"
       
    91 ; delete(clat@units)
       
    92 ; printVarSummary(clat)
       
    93 
       
    94  bo = new((/12,nlat,nlon/),float)
       
    95 ;bo = new((/1,nlat,nlon/),float)
       
    96  do m = 0,11
       
    97     do j=0,nlat-1
       
    98        if (j.eq.0 .or. j.eq.nlat-1) then
       
    99           if (j.eq.0) then
       
   100              LATS = -90.          
       
   101              LATN = lat(j)+0.5*(lat(j+1)-lat(j))
       
   102           end if
       
   103           if (j.eq.nlat-1) then
       
   104              LATS = lat(j)-0.5*(lat(j)-lat(j-1))
       
   105              LATN = 90.                  
       
   106           end if
       
   107        else
       
   108           LATS = lat(j)-0.5*(lat(j)-lat(j-1))
       
   109           LATN = lat(j)+0.5*(lat(j+1)-lat(j))
       
   110        end if
       
   111  
       
   112 ;      CLAT = clat({LATS:LATN})      ; do once for *slight* efficiency
       
   113 ;      TEMP = bi(:,{LATS:LATN},:)      ; 2D [lat,lon]
       
   114  
       
   115       do i=0,nlon-1
       
   116        if (i.eq.0 .or. i.eq.nlon-1) then
       
   117           if (i.eq.0) then
       
   118              LONL = 0.          
       
   119              LONR = lon(i)+0.5*(lon(i+1)-lon(i))
       
   120           end if
       
   121           if (i.eq.nlon-1) then
       
   122              LONL = lon(i)-0.5*(lon(i)-lon(i-1))
       
   123              LONR = 360.                 
       
   124           end if
       
   125        else
       
   126           LONL = lon(i)-0.5*(lon(i)-lon(i-1))
       
   127           LONR = lon(i)+0.5*(lon(i+1)-lon(i))
       
   128        end if
       
   129 
       
   130 ;print (LATS)
       
   131 ;print (LATN)
       
   132 ;print (LONL)
       
   133 ;print (LONR)
       
   134 ;        bo(:,j,i) = avg(bi(:,{LATS:LATN},{LONL:LONR}))
       
   135          bo(m,j,i) = avg(bi(m,{LATS:LATN},{LONL:LONR}))  
       
   136 ;        bo(:,j,i) = wgt_areaave(TEMP(:,{LONL:LONR}), CLAT, 1.0, 0)
       
   137       end do
       
   138  
       
   139 ;     delete(CLAT)
       
   140 ;     delete(TEMP) 
       
   141     end do
       
   142   end do
       
   143 
       
   144   bo!0   = "time"
       
   145   bo!1   = "lat"
       
   146   bo!2   = "lon"
       
   147   bo&time= bi&time
       
   148 ; bo&time= 1
       
   149   bo&lat = lat
       
   150   bo&lon = lon
       
   151 ; bo@units      = bi@units
       
   152 ; bo@long_name  = bi@long_name
       
   153   bo@units      = "none"
       
   154   bo@long_name  = "Leaf Area Index"
       
   155   bo@_FillValue = bi@_FillValue
       
   156 
       
   157   c->LAI  = bo
       
   158 end