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/31.contour_ob_grow.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 ;*************************************************

     2 ; ce_1.ncl

     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

     9 ;************************************************

    10 ; read in observed data

    11 ;************************************************

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

    13   fili  = "LAI_2000-2005_ensemble_T42.nc"

    14   f     = addfile(diri+fili,"r")

    15 

    16   z     = f->LAI

    17   y     = z(0,:,:)

    18   y@long_name = "Days of Growing Season"

    19  

    20   dsizes_z = dimsizes(z)

    21   ntime    = dsizes_z(0)

    22   nlat     = dsizes_z(1)

    23   nlon     = dsizes_z(2)

    24 

    25  day_of_data = (/31,28,31,30,31,30,31,31,30,31,30,31/)

    26 

    27   do j = 0,nlat-1

    28   do i = 0,nlon-1

    29      nday = 0.

    30      do k = 0,ntime-1

    31         if (.not. ismissing(z(k,j,i)) .and. z(k,j,i) .gt. 1.0) then

    32            nday = nday + day_of_data(k)

    33         end if

    34      end do

    35      y(j,i) = nday

    36   end do

    37   end do

    38 

    39   print (min(y)+"/"+max(y))

    40 

    41 ;************************************************

    42 ; read in model data

    43 ;************************************************

    44   diri2  = "/fis/cgd/cseg/people/jeff/clamp_data/model/"

    45   fili2 = "i01.03cn_1545-1569_ANN_climo.nc"

    46   g     = addfile(diri2+fili2,"r")

    47   x     = g->TLAI

    48 

    49   delta = 0.000001

    50   x0    = x(0,:,:)

    51   y = where(ismissing(y).and.(ismissing(x0).or.(x0.lt.delta)),y@_FillValue,y)

    52 ;************************************************

    53 ; create default plot

    54 ;************************************************

    55 

    56   wks = gsn_open_wks("ps","xy")         ; open a ps file

    57   gsn_define_colormap(wks,"gui_default")     ; choose colormap

    58 

    59   res                     = True             ; Use plot options

    60   res@cnFillOn            = True             ; Turn on color fill

    61   res@gsnSpreadColors      = True            ; use full colormap

    62 ; res@cnFillMode          = "RasterFill"     ; Turn on raster color

    63 ; res@lbLabelAutoStride   = True

    64   res@cnLinesOn           = False            ; Turn off contourn lines

    65   res@mpFillOn            = False            ; Turn off map fill

    66   res@tiMainString        = "MODIS MOD 15A2 2000-2005"

    67 

    68   res@gsnSpreadColors     = True             ; use full colormap

    69   res@cnLevelSelectionMode = "ManualLevels"  ; Manual contour invtervals

    70   res@cnMinLevelValF       = 0.              ; Min level

    71   res@cnMaxLevelValF       = 390.            ; Max level

    72   res@cnLevelSpacingF      = 30.             ; interval

    73 

    74 ; pres                            = True     ; panel plot mods desired

    75 ; pres@gsnMaximize                = True     ; fill the page

    76 

    77   plot = gsn_csm_contour_map_ce(wks,y,res)   ; for observed

    78 

    79 ; gsn_panel(wks,plot,(/1,1/),pres)           ; create panel plot

    80   system("convert xy.ps xy.png")

    81 end
clamp