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