diff -r 000000000000 -r 0c6405ab2ff4 lai/33.contour_diff_mean.ncl
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/lai/33.contour_diff_mean.ncl	Mon Jan 26 22:08:20 2009 -0500
@@ -0,0 +1,114 @@
+;*************************************************
+; 3 plots, model - ob
+;************************************************
+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
+;************************************************
+  diri1  = "/fis/cgd/cseg/people/jeff/clamp_data/lai/"
+  fili1  = "LAI_2000-2005_ensemble_T42.nc"
+  f     = addfile(diri1+fili1,"r")
+
+  z     = f->LAI
+  y     = z(0,:,:)
+  y@long_name = "Leaf Area Index Mean"
+  s     = z(:,0,0)
+ 
+  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) = avg(s)
+  end do
+  end do
+
+; printVarSummary(y)
+  print (min(y)+"/"+max(y))
+
+  delete (s)
+  delete (z)
+;************************************************
+; read in data: model       
+;************************************************
+ diri2  = "/fis/cgd/cseg/people/jeff/clamp_data/model/"
+ fili3  = "i01.03cn_1545-1569_MONS_climo.nc"
+;fili3  = "i01.04casa_1605-1629_MONS_climo.nc"
+ data_file_model = addfile(diri2+fili3,"r")
+
+  z     = data_file_model->TLAI
+  x     = z(0,:,:)
+  x@long_name = "Leaf Area Index Mean"
+  s     = z(:,0,0)  
+ 
+  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) 
+     x(j,i) = avg(s)
+  end do
+  end do
+
+  print (min(x)+"/"+max(x))
+  d = x              
+  d = x - y
+  print (min(d)+"/"+max(d))
+
+  delete (z)
+  delete (s)
+
+  delta = 0.000001
+  y = where(ismissing(y).and.(ismissing(x).or.(x.lt.delta)),0.,y)
+;************************************************
+; create default plot
+;************************************************
+  wks = gsn_open_wks("ps","xy")              ; open a ps file
+  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@gsnSpreadColors      = True            ; use full colormap
+  res@cnLevelSelectionMode = "ManualLevels"  ; Manual contour invtervals
+  res@cnMinLevelValF       = 0.              ; Min level
+  res@cnMaxLevelValF       = 6.             ; Max level
+  res@cnLevelSpacingF      = 0.5             ; interval
+
+  pres                     = True            ; panel plot mods desired
+  pres@gsnMaximize         = True            ; fill the page
+
+  plot=new(3,graphic)                        ; create graphic array
+
+  res@tiMainString        = "MODIS MOD 15A2 2000-2005"
+  plot(0) = gsn_csm_contour_map_ce(wks,y,res)   
+
+  res@tiMainString        = "Model i01.03cn"
+; res@tiMainString        = "Model i01.04casa"
+  plot(1) = gsn_csm_contour_map_ce(wks,x,res)   
+ 
+  res@tiMainString        = "(Model i01.03cn) - (Observed)"
+; res@tiMainString        = "(Model i01.04casa) - (Observed)"    
+  res@cnMinLevelValF       = -2.              ; Min level
+  res@cnMaxLevelValF       =  2.             ; Max level
+  res@cnLevelSpacingF      =  0.4              ; interval
+  plot(2) = gsn_csm_contour_map_ce(wks,d,res)   ; for observed
+
+  gsn_panel(wks,plot,(/3,1/),pres)           ; create panel plot
+  system("convert xy.ps xy.png")
+end
\ No newline at end of file