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.
forrest@0
     1
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"  
forrest@0
     2
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"  
forrest@0
     3
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"  
forrest@0
     4
forrest@0
     5
begin 
forrest@0
     6
;----------------------------------------------------------
forrest@0
     7
  year = "2004"
forrest@0
     8
  
forrest@0
     9
  nlat  = 3600
forrest@0
    10
  mlon  = 7200
forrest@0
    11
  month_in_year = 12
forrest@0
    12
;----------------------------------------------------------
forrest@0
    13
  diri  = "/fis/cgd/cseg/people/jeff/clamp_data/"
forrest@0
    14
  diro  = "/fis/cgd/cseg/people/jeff/clamp_data/"
forrest@0
    15
  filo  = "MOD15A2_LAI_" + year + "_monthly.nc"
forrest@0
    16
  c = addfile(diro+filo,"c")
forrest@0
    17
  filedimdef(c,"time",-1,True)
forrest@0
    18
forrest@0
    19
  day = (/"001","009","017","025","033","041","049","057","065","073", \
forrest@0
    20
          "081","089","097","105","113","121","129","137","145","153", \
forrest@0
    21
          "161","169","177","185","193","201","209","217","225","233", \
forrest@0
    22
          "241","249","257","265","273","281","289","297","305","313", \
forrest@0
    23
          "321","329","337","345","353","361"/)
forrest@0
    24
forrest@0
    25
  file_index_L = (/0,3,7 ,11,15,18,22,26,30,34,38,41/)
forrest@0
    26
  file_index_R = (/3,7,11,14,18,22,26,30,34,37,41,45/)
forrest@0
    27
  scale_L      = (/8,1,5,6,8,1,3,4,5,7,8,2/)
forrest@0
    28
  scale_R      = (/7,3,2,8,7,5,4,3,1,8,6,5/)
forrest@0
    29
forrest@0
    30
;  dayI = stringtointeger(day)
forrest@0
    31
;  print (dayI)
forrest@0
    32
forrest@0
    33
  x = new((/12,nlat,mlon/),float)
forrest@0
    34
  x = 0.
forrest@0
    35
  x@_FillValue =  1.e+36
forrest@0
    36
forrest@0
    37
 do m = 0,month_in_year-1
forrest@0
    38
    nday = 0                                   ; number of day in a month
forrest@0
    39
 do n = file_index_L(m),file_index_R(m)
forrest@0
    40
    fili  = "MOD15A2_GEO_" + year + day(n) + ".hdf"
forrest@0
    41
    print (fili)
forrest@0
    42
    a = addfile(diri+fili,"r")
forrest@0
    43
    y = a->Lai_0_05deg
forrest@0
    44
forrest@0
    45
    y@_FillValue = inttobyte(253)              ; special missing value? 
forrest@0
    46
    y@_FillValue = inttobyte(255)              ; missing value 
forrest@0
    47
forrest@0
    48
    z = byte2flt(y)
forrest@0
    49
          
forrest@0
    50
    z@_FillValue =  1.e+36
forrest@0
    51
forrest@0
    52
;   weighted by day, data is 8-day average                 
forrest@0
    53
    scale = 8              
forrest@0
    54
    if (n .eq. file_index_L(m)) then
forrest@0
    55
       scale = scale_L(m)
forrest@0
    56
    end if
forrest@0
    57
    if (n .eq. file_index_R(m)) then
forrest@0
    58
       scale = scale_R(m)
forrest@0
    59
    end if
forrest@0
    60
forrest@0
    61
    nday = nday + scale
forrest@0
    62
forrest@0
    63
    x(m,:,:) = where(ismissing(z),1.e+36,x(m,:,:)+z(:,:)*scale)
forrest@0
    64
forrest@0
    65
    delete (a)
forrest@0
    66
    delete (y)
forrest@0
    67
    delete (z)
forrest@0
    68
 end do
forrest@0
    69
    x(m,:,:) = where(ismissing(x(m,:,:)),1.e+36,x(m,:,:)/nday)
forrest@0
    70
    print (nday)
forrest@0
    71
    print (min(x) + "/" + max(x))
forrest@0
    72
 end do
forrest@0
    73
forrest@0
    74
  lat  = latGlobeFo(nlat, "lat", "latitude", "degrees_north")
forrest@0
    75
  lat  = (/ lat(::-1) /)                     ; make N->S
forrest@0
    76
  lon  = lonGlobeFo(mlon, "lon", "longitude", "degrees_east")
forrest@0
    77
  lon  = (/ lon - 180. /) ; subtract 180 from all values 
forrest@0
    78
  lon&lon = lon           ; update coordinates
forrest@0
    79
forrest@0
    80
  x!0  = "time"
forrest@0
    81
  x!1  = "lat"
forrest@0
    82
  x!2  = "lon"
forrest@0
    83
  x&time= ispan(1,month_in_year,1)
forrest@0
    84
  x&lat=  lat
forrest@0
    85
  x&lon=  lon
forrest@0
    86
  x@units      = "none"
forrest@0
    87
  x@long_name  = "Leaf Area Index"
forrest@0
    88
forrest@0
    89
  c->LAI  = x 
forrest@0
    90
end
forrest@0
    91