Mercurial Mercurial > repos > public_software > clamp / file revision
summary | shortlog | changelog | graph | tags | bookmarks | branches | files | changeset | file | latest | revisions | annotate | diff | comparison | raw | help
lai/02.read_lai_avg.ncl
author Forrest Hoffman <forrest@climatemodeling.org>
Thu, 26 Mar 2009 14:02:21 -0400
changeset 1 4be95183fbcd
permissions -rw-r--r--
Modifications to scoring and graphics production for the final version of code for the C-LAMP paper in GCB. 
     1 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"  

     2 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"  

     3 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"  

     4 

     5 begin 

     6 ;----------------------------------------------------------

     7   year = "2004"

     8   

     9   nlat  = 3600

    10   mlon  = 7200

    11   month_in_year = 12

    12 ;----------------------------------------------------------

    13   diri  = "/fis/cgd/cseg/people/jeff/clamp_data/"

    14   diro  = "/fis/cgd/cseg/people/jeff/clamp_data/"

    15   filo  = "MOD15A2_LAI_" + year + "_monthly.nc"

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

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

    18 

    19   day = (/"001","009","017","025","033","041","049","057","065","073", \

    20           "081","089","097","105","113","121","129","137","145","153", \

    21           "161","169","177","185","193","201","209","217","225","233", \

    22           "241","249","257","265","273","281","289","297","305","313", \

    23           "321","329","337","345","353","361"/)

    24 

    25   file_index_L = (/0,3,7 ,11,15,18,22,26,30,34,38,41/)

    26   file_index_R = (/3,7,11,14,18,22,26,30,34,37,41,45/)

    27   scale_L      = (/8,1,5,6,8,1,3,4,5,7,8,2/)

    28   scale_R      = (/7,3,2,8,7,5,4,3,1,8,6,5/)

    29 

    30 ;  dayI = stringtointeger(day)

    31 ;  print (dayI)

    32 

    33   x = new((/12,nlat,mlon/),float)

    34   x = 0.

    35   x@_FillValue =  1.e+36

    36 

    37  do m = 0,month_in_year-1

    38     nday = 0                                   ; number of day in a month

    39  do n = file_index_L(m),file_index_R(m)

    40     fili  = "MOD15A2_GEO_" + year + day(n) + ".hdf"

    41     print (fili)

    42     a = addfile(diri+fili,"r")

    43     y = a->Lai_0_05deg

    44 

    45     y@_FillValue = inttobyte(253)              ; special missing value? 

    46     y@_FillValue = inttobyte(255)              ; missing value 

    47 

    48     z = byte2flt(y)

    49           

    50     z@_FillValue =  1.e+36

    51 

    52 ;   weighted by day, data is 8-day average                 

    53     scale = 8              

    54     if (n .eq. file_index_L(m)) then

    55        scale = scale_L(m)

    56     end if

    57     if (n .eq. file_index_R(m)) then

    58        scale = scale_R(m)

    59     end if

    60 

    61     nday = nday + scale

    62 

    63     x(m,:,:) = where(ismissing(z),1.e+36,x(m,:,:)+z(:,:)*scale)

    64 

    65     delete (a)

    66     delete (y)

    67     delete (z)

    68  end do

    69     x(m,:,:) = where(ismissing(x(m,:,:)),1.e+36,x(m,:,:)/nday)

    70     print (nday)

    71     print (min(x) + "/" + max(x))

    72  end do

    73 

    74   lat  = latGlobeFo(nlat, "lat", "latitude", "degrees_north")

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

    76   lon  = lonGlobeFo(mlon, "lon", "longitude", "degrees_east")

    77   lon  = (/ lon - 180. /) ; subtract 180 from all values 

    78   lon&lon = lon           ; update coordinates

    79 

    80   x!0  = "time"

    81   x!1  = "lat"

    82   x!2  = "lon"

    83   x&time= ispan(1,month_in_year,1)

    84   x&lat=  lat

    85   x&lon=  lon

    86   x@units      = "none"

    87   x@long_name  = "Leaf Area Index"

    88 

    89   c->LAI  = x 

    90 end

    91
clamp