1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/lai/03.land_class_to_T42.ncl Mon Jan 26 22:08:20 2009 -0500
1.3 @@ -0,0 +1,141 @@
1.4 +; ***********************************************
1.5 +; interpolate into model grids (T42)
1.6 +; output: lat: S->N lon: 0->360
1.7 +; input : lat: N->S lon: -180->180
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 + nclass = 20
1.16 +
1.17 +;************************************************
1.18 +; output data
1.19 +;************************************************
1.20 + diro = "/fis/cgd/cseg/people/jeff/clamp_data/lai/"
1.21 +; filo = "land_class_T42.nc"
1.22 + filo = "land_class_T42_new.nc"
1.23 + c = addfile(diro+filo,"c")
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 = "land_class_0.25deg.nc"
1.30 + fili = "land_class_0.25deg_new.nc"
1.31 + g = addfile (diri+fili,"r")
1.32 + bi = g->LAND_CLASS
1.33 + xi = g->lon
1.34 + yi = g->lat
1.35 +
1.36 +;************************************************
1.37 +; change into 0-360E, 90S-90N
1.38 +;************************************************
1.39 +
1.40 + yi = (/ yi(::-1) /)
1.41 + bi = (/ bi(::-1,:) /)
1.42 + printVarSummary(bi)
1.43 +
1.44 + b2 = bi
1.45 + x2 = xi
1.46 +
1.47 + nx = dimsizes(xi)
1.48 + nx2 = nx/2
1.49 +
1.50 + do i= 0,nx-1
1.51 + if (i .lt. nx2) then
1.52 + p = i + nx2
1.53 + xi(p) = x2(i) + 360.
1.54 + else
1.55 + p = i - nx2
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 +; read in model data
1.70 +;************************************************
1.71 + diri2 = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
1.72 + fili2 = "i01.03cn_1545-1569_ANN_climo.nc"
1.73 + f = addfile (diri2+fili2,"r")
1.74 +
1.75 + lon = f->lon
1.76 + lat = f->lat
1.77 + nlon = dimsizes(lon)
1.78 + nlat = dimsizes(lat)
1.79 +
1.80 +; print (xi)
1.81 +; print (yi)
1.82 +; print (xo)
1.83 +; print (yo)
1.84 +
1.85 + bo = new((/nlat,nlon/),integer)
1.86 + count = new((/nclass/),integer)
1.87 +
1.88 + do j=0,nlat-1
1.89 + if (j.eq.0 .or. j.eq.nlat-1) then
1.90 + if (j.eq.0) then
1.91 + LATS = -90.
1.92 + LATN = lat(j)+0.5*(lat(j+1)-lat(j))
1.93 + end if
1.94 + if (j.eq.nlat-1) then
1.95 + LATS = lat(j)-0.5*(lat(j)-lat(j-1))
1.96 + LATN = 90.
1.97 + end if
1.98 + else
1.99 + LATS = lat(j)-0.5*(lat(j)-lat(j-1))
1.100 + LATN = lat(j)+0.5*(lat(j+1)-lat(j))
1.101 + end if
1.102 +
1.103 + do i=0,nlon-1
1.104 + if (i.eq.0 .or. i.eq.nlon-1) then
1.105 + if (i.eq.0) then
1.106 + LONL = 0.
1.107 + LONR = lon(i)+0.5*(lon(i+1)-lon(i))
1.108 + end if
1.109 + if (i.eq.nlon-1) then
1.110 + LONL = lon(i)-0.5*(lon(i)-lon(i-1))
1.111 + LONR = 360.
1.112 + end if
1.113 + else
1.114 + LONL = lon(i)-0.5*(lon(i)-lon(i-1))
1.115 + LONR = lon(i)+0.5*(lon(i+1)-lon(i))
1.116 + end if
1.117 +
1.118 +;print (LATS)
1.119 +;print (LATN)
1.120 +;print (LONL)
1.121 +;print (LONR)
1.122 +
1.123 + count = 0
1.124 +
1.125 + do n = 0,nclass-1
1.126 + count(n) = num (bi({LATS:LATN},{LONL:LONR}).eq.n)
1.127 + end do
1.128 +
1.129 + bo(j,i) = maxind(count)
1.130 + end do
1.131 + end do
1.132 +
1.133 + bo!0 = "lat"
1.134 + bo!1 = "lon"
1.135 + bo&lat = lat
1.136 + bo&lon = lon
1.137 + bo@units = bi@units
1.138 + bo@long_name = bi@long_name
1.139 + bo@_FillValue = bi@_FillValue
1.140 +
1.141 + print (min(bo) + "/" + max(bo))
1.142 +
1.143 + c->LAND_CLASS = bo
1.144 +end
1.145 \ No newline at end of file