clamp: lai/03.lai_ob_to_0.25deg.ncl@4be95183fbcd

     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

    10   year = 2000

    12 ;************************************************

    13 ; output data

    14 ;************************************************

    15   diro  = "/fis/cgd/cseg/people/jeff/clamp_data/lai/"

    16 ; filo  = "LAI_" + year + "_monthly_0.25deg.nc"

    17   filo  = "LAI_2000-2005_mean_0.25deg.nc"

    18 ; filo  = "LAI_2000-2005_ensemble_0.25deg.nc"

    19   c = addfile(diro+filo,"c")

    20   filedimdef(c,"time",-1,True)

    22 ;************************************************

    23 ; read in observed data

    24 ;************************************************

    25   diri  = "/fis/cgd/cseg/people/jeff/clamp_data/lai/"

    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

    34 ;************************************************

    35 ; change into 0-360E, 90S-90N

    36 ; Observed NPP*scale_factor

    37 ;************************************************

    39  yi    = (/ yi(::-1) /)

    40  bi    = (/ bi(:,::-1,:) /)

    41  printVarSummary(bi)

    43  b2    = bi

    44  x2    = xi

    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

    58  bi&lat =  yi

    59  bi&lon =  xi

    61 ;print (xi)

    62 ;print (yi)

    63 ;exit

    65 ;************************************************

    66 ; create 0.25deg lat and lon

    67 ;************************************************

    68   nlon   = 360*4

    69   nlat   = 180*4

    71   lat  = latGlobeFo(nlat, "lat", "latitude", "degrees_north") ; S->N

    72 ; lat  = (/ lat(::-1) /)                                      ; N->S

    73   lat&lat = lat

    74   lon  = lonGlobeFo(nlon, "lon", "longitude", "degrees_east") ; 0->360

    75 ; lon  = (/ lon - 180. /) ; subtract 180 from all values      ; 180W-180E

    76   lon&lon = lon           ; update coordinates

    78 ; print (lon)

    79 ; print (lat)

    81 ; rad   = 4.*atan(1.)/180.

    82 ; clat  = lat

    83 ; clat  = cos(lat*rad)

    84 ; clat@long_name = "cos(latitude)"

    85 ; delete(clat@units)

    86 ; printVarSummary(clat)

    88 ;bo = new((/12,nlat,nlon/),float)

    89  bo = new((/1,nlat,nlon/),float)

    91 ;do m = 0,11

    92  do m = 0,0

    93     do j=0,nlat-1

    94        if (j.eq.0 .or. j.eq.nlat-1) then

    95           if (j.eq.0) then

    96              LATS = -90.

    97              LATN = lat(j)+0.5*(lat(j+1)-lat(j))

    98           end if

    99           if (j.eq.nlat-1) then

   100              LATS = lat(j)-0.5*(lat(j)-lat(j-1))

   101              LATN = 90.

   102           end if

   103        else

   104           LATS = lat(j)-0.5*(lat(j)-lat(j-1))

   105           LATN = lat(j)+0.5*(lat(j+1)-lat(j))

   106        end if

   108 ;      CLAT = clat({LATS:LATN})      ; do once for *slight* efficiency

   109 ;      TEMP = bi(:,{LATS:LATN},:)      ; 2D [lat,lon]

   111       do i=0,nlon-1

   112        if (i.eq.0 .or. i.eq.nlon-1) then

   113           if (i.eq.0) then

   114              LONL = 0.

   115              LONR = lon(i)+0.5*(lon(i+1)-lon(i))

   116           end if

   117           if (i.eq.nlon-1) then

   118              LONL = lon(i)-0.5*(lon(i)-lon(i-1))

   119              LONR = 360.

   120           end if

   121        else

   122           LONL = lon(i)-0.5*(lon(i)-lon(i-1))

   123           LONR = lon(i)+0.5*(lon(i+1)-lon(i))

   124        end if

   126 ;print (LATS)

   127 ;print (LATN)

   128 ;print (LONL)

   129 ;print (LONR)

   131          bo(m,j,i) = avg(bi(m,{LATS:LATN},{LONL:LONR}))

   132 ;        bo(m,j,i) = wgt_areaave(TEMP(m,{LONL:LONR}), CLAT, 1.0, 0)

   133       end do

   135 ;     delete(CLAT)

   136 ;     delete(TEMP)

   137     end do

   138   end do

   140   bo!0   = "time"

   141   bo!1   = "lat"

   142   bo!2   = "lon"

   143   bo&time= bi&time

   144 ; bo&time= 1

   145   bo&lat = lat

   146   bo&lon = lon

   147 ; bo@units      = bi@units

   148 ; bo@long_name  = bi@long_name

   149   bo@units      = "none"

   150   bo@long_name  = "Leaf Area Index"

   151   bo@_FillValue = bi@_FillValue

   153   c->LAI  = bo

   154 end

author	Forrest Hoffman <forrest@climatemodeling.org>
	Thu, 26 Mar 2009 14:02:21 -0400
changeset 1	4be95183fbcd
permissions	-rw-r--r--