diff -r 000000000000 -r 0c6405ab2ff4 lai/31.contour_ob_max.ncl
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/lai/31.contour_ob_max.ncl	Mon Jan 26 22:08:20 2009 -0500
@@ -0,0 +1,98 @@
+;*************************************************
+; ce_1.ncl
+;************************************************
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"  
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"  
+;************************************************
+begin
+;************************************************
+; read in observed data
+;************************************************
+  diri  = "/fis/cgd/cseg/people/jeff/clamp_data/"
+  fili  = "LAI_2000-2005_ensemble_T42.nc"
+; fili  = "LAI_2000_monthly_T42.nc"
+  f     = addfile(diri+fili,"r")
+
+  z     = f->LAI
+  y     = z(0,:,:)
+  y@long_name = "Leaf Area Index Max"
+  p     = z(0,:,:)
+  p@long_name = "Leaf Area Index Max Month"
+  s     = z(:,0,0)
+; printVarSummary(s)
+ 
+  dsizes_z = dimsizes(z)
+  ntime    = dsizes_z(0)
+  nlat     = dsizes_z(1)
+  nlon     = dsizes_z(2)
+  
+  do j = 0,nlat-1
+  do i = 0,nlon-1
+     s      = z(:,j,i) 
+     y(j,i) = max(s)
+     p(j,i) = maxind(s)+ 1
+
+;    iMax   = ind(s.eq.y(j,i)) + 1      ; may be multiple max
+;    if (.not.ismissing(iMax(0))) then 
+;        p(j,i) =  iMax(0)          ; no meta data
+;    end if
+;    delete(iMax)
+
+  end do
+  end do
+
+; printVarSummary(p)
+  print (min(y)+"/"+max(y))
+  print (min(p)+"/"+max(p))
+
+;************************************************
+; read in model data
+;************************************************
+  fili2 = "i01.03cn_1545-1569_ANN_climo.nc"
+  g     = addfile(diri+fili2,"r")
+  x     = g->TLAI
+
+  delta = 0.00000000001
+  x0    = x(0,:,:)
+  y = where(ismissing(y).and.(ismissing(x0).or.(x0.lt.delta)),0.,y)
+;************************************************
+; create default plot
+;************************************************
+  
+; setvalues NhlGetWorkspaceObjectId()
+;   "wsMaximumSize" : 199999999
+; end setvalues
+
+  wks = gsn_open_wks("ps","xy")         ; open a ps file
+; gsn_define_colormap(wks,"wgne15")          ; choose colormap
+  gsn_define_colormap(wks,"gui_default")     ; choose colormap
+
+  res                     = True             ; Use plot options
+  res@cnFillOn            = True             ; Turn on color fill
+  res@gsnSpreadColors      = True            ; use full colormap
+; res@cnFillMode          = "RasterFill"     ; Turn on raster color
+; res@lbLabelAutoStride   = True
+  res@cnLinesOn           = False            ; Turn off contourn lines
+  res@mpFillOn            = False            ; Turn off map fill
+  res@tiMainString        = "MODIS MOD 15A2 2000-2005"
+
+  res@gsnSpreadColors     = True             ; use full colormap
+  res@cnLevelSelectionMode = "ManualLevels"  ; Manual contour invtervals
+  res@cnMinLevelValF       = 0.              ; Min level
+  res@cnMaxLevelValF       = 10.             ; Max level
+  res@cnLevelSpacingF      = 1.              ; interval
+
+; pres                            = True     ; panel plot mods desired
+; pres@gsnMaximize                = True     ; fill the page
+
+  plot = gsn_csm_contour_map_ce(wks,y,res)   ; for observed
+
+  res@cnMinLevelValF       = 1.              ; Min level
+  res@cnMaxLevelValF       = 12.             ; Max level
+  res@cnLevelSpacingF      = 1.              ; interval
+  plot = gsn_csm_contour_map_ce(wks,p,res)   ; for observed
+
+; gsn_panel(wks,plot,(/1,1/),pres)           ; create panel plot
+  system("convert xy.ps xy.png")
+end
\ No newline at end of file