lai/03.lai_ob_to_0.25deg.ncl.0
changeset 0 0c6405ab2ff4
     1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/lai/03.lai_ob_to_0.25deg.ncl.0	Mon Jan 26 22:08:20 2009 -0500
     1.3 @@ -0,0 +1,157 @@
     1.4 +; ***********************************************
     1.5 +; interpolate into 0.25degree
     1.6 +; no change in lat and lon, to be in consistent with 
     1.7 +; land class data
     1.8 +; ***********************************************
     1.9 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
    1.10 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
    1.11 +load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
    1.12 +;************************************************
    1.13 +begin
    1.14 +
    1.15 +  year = 2003
    1.16 +
    1.17 +;************************************************
    1.18 +; output data
    1.19 +;************************************************
    1.20 +  diro  = "/fis/cgd/cseg/people/jeff/clamp_data/lai/"
    1.21 +  filo  = "LAI_2000-2005_mean_0.25.nc"
    1.22 +  c = addfile(diro+filo,"c")
    1.23 +  filedimdef(c,"time",-1,True)
    1.24 +
    1.25 +;************************************************
    1.26 +; read in observed data
    1.27 +;************************************************
    1.28 +  diri  = "/fis/cgd/cseg/people/jeff/clamp_data/lai/"
    1.29 +; fili  = "LAI_2000-2005_mean_0.05deg.nc"
    1.30 +  fili  = "LAI_2000-2005_ensemble_0.05deg.nc"
    1.31 +  g     = addfile (diri+fili,"r")
    1.32 +  bi    = g->LAI   
    1.33 +  xi    = g->lon 
    1.34 +  yi    = g->lat
    1.35 +
    1.36 +;************************************************
    1.37 +; change into 0-360E, 90S-90N
    1.38 +; Observed NPP*scale_factor
    1.39 +;************************************************
    1.40 + 
    1.41 + yi    = (/ yi(::-1) /)
    1.42 + bi    = (/ bi(:,::-1,:) /)
    1.43 + printVarSummary(bi)
    1.44 +
    1.45 + b2    = bi
    1.46 + x2    = xi   
    1.47 + 
    1.48 + nx = dimsizes(xi)
    1.49 + do i= 0,nx-1
    1.50 +    if (i .lt. 3600) then
    1.51 +       p = i + 3600
    1.52 +       xi(p) = x2(i) + 360.      
    1.53 +    else
    1.54 +       p = i - 3600
    1.55 +       xi(p) = x2(i)
    1.56 +    end if
    1.57 +    bi(:,:,p)= b2(:,:,i) 
    1.58 + end do
    1.59 +
    1.60 + bi&lat =  yi
    1.61 + bi&lon =  xi
    1.62 +
    1.63 +;print (xi)
    1.64 +;print (yi)
    1.65 +;exit
    1.66 +
    1.67 +;************************************************
    1.68 +; read in model data
    1.69 +;************************************************
    1.70 +;fili2  = "i01.03cn_1545-1569_ANN_climo.nc"
    1.71 + fili2  = ""
    1.72 + f     = addfile (diri+fili2,"r")
    1.73 +
    1.74 + lon    = f->lon     
    1.75 + lat    = f->lat
    1.76 + nlon   = dimsizes(lon)
    1.77 + nlat   = dimsizes(lat)      
    1.78 +
    1.79 +; print (xi)
    1.80 +; print (yi)
    1.81 +; print (xo)
    1.82 +; print (yo)
    1.83 +
    1.84 +;(A)
    1.85 +;bo = linint2_Wrap(xi,yi,bi,True,xo,yo,0)
    1.86 +
    1.87 +;(B)
    1.88 +
    1.89 +; rad   = 4.*atan(1.)/180.
    1.90 +; clat  = lat
    1.91 +; clat  = cos(lat*rad)
    1.92 +; clat@long_name = "cos(latitude)"
    1.93 +; delete(clat@units)
    1.94 +; printVarSummary(clat)
    1.95 +
    1.96 + bo = new((/12,nlat,nlon/),float)
    1.97 +;bo = new((/1,nlat,nlon/),float)
    1.98 + do m = 0,11
    1.99 +    do j=0,nlat-1
   1.100 +       if (j.eq.0 .or. j.eq.nlat-1) then
   1.101 +          if (j.eq.0) then
   1.102 +             LATS = -90.          
   1.103 +             LATN = lat(j)+0.5*(lat(j+1)-lat(j))
   1.104 +          end if
   1.105 +          if (j.eq.nlat-1) then
   1.106 +             LATS = lat(j)-0.5*(lat(j)-lat(j-1))
   1.107 +             LATN = 90.                  
   1.108 +          end if
   1.109 +       else
   1.110 +          LATS = lat(j)-0.5*(lat(j)-lat(j-1))
   1.111 +          LATN = lat(j)+0.5*(lat(j+1)-lat(j))
   1.112 +       end if
   1.113 + 
   1.114 +;      CLAT = clat({LATS:LATN})      ; do once for *slight* efficiency
   1.115 +;      TEMP = bi(:,{LATS:LATN},:)      ; 2D [lat,lon]
   1.116 + 
   1.117 +      do i=0,nlon-1
   1.118 +       if (i.eq.0 .or. i.eq.nlon-1) then
   1.119 +          if (i.eq.0) then
   1.120 +             LONL = 0.          
   1.121 +             LONR = lon(i)+0.5*(lon(i+1)-lon(i))
   1.122 +          end if
   1.123 +          if (i.eq.nlon-1) then
   1.124 +             LONL = lon(i)-0.5*(lon(i)-lon(i-1))
   1.125 +             LONR = 360.                 
   1.126 +          end if
   1.127 +       else
   1.128 +          LONL = lon(i)-0.5*(lon(i)-lon(i-1))
   1.129 +          LONR = lon(i)+0.5*(lon(i+1)-lon(i))
   1.130 +       end if
   1.131 +
   1.132 +;print (LATS)
   1.133 +;print (LATN)
   1.134 +;print (LONL)
   1.135 +;print (LONR)
   1.136 +;        bo(:,j,i) = avg(bi(:,{LATS:LATN},{LONL:LONR}))
   1.137 +         bo(m,j,i) = avg(bi(m,{LATS:LATN},{LONL:LONR}))  
   1.138 +;        bo(:,j,i) = wgt_areaave(TEMP(:,{LONL:LONR}), CLAT, 1.0, 0)
   1.139 +      end do
   1.140 + 
   1.141 +;     delete(CLAT)
   1.142 +;     delete(TEMP) 
   1.143 +    end do
   1.144 +  end do
   1.145 +
   1.146 +  bo!0   = "time"
   1.147 +  bo!1   = "lat"
   1.148 +  bo!2   = "lon"
   1.149 +  bo&time= bi&time
   1.150 +; bo&time= 1
   1.151 +  bo&lat = lat
   1.152 +  bo&lon = lon
   1.153 +; bo@units      = bi@units
   1.154 +; bo@long_name  = bi@long_name
   1.155 +  bo@units      = "none"
   1.156 +  bo@long_name  = "Leaf Area Index"
   1.157 +  bo@_FillValue = bi@_FillValue
   1.158 +
   1.159 +  c->LAI  = bo
   1.160 +end
   1.161 \ No newline at end of file